AFIT/EN /EN  G  /  97D-01 


ANALYTIC  TRANSFER  FUNCTION  OF  THE 
FORWARD  PROPAGATION  OF  DIFFUSE  PHOTON 
DENSITY  WAVES  IN  TURBID  MEDIA  WITH 
AN  EMBEDDED  SPHERICAL  INHOMOGENEITY 

THESIS 

Deborah  L.  Lasocki 
Captain,  USAF 

AFIT/EN/ENG/97D-01 


19980128  111 


Approved  for  public  release;  distribution  unlimited 


\jma  QUALITY-  BJgPZOTED  3 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  Department  of  Defense  or  the  United  States  Government. 


AFIT/EN/ENG/97D-01 


ANALYTIC  TRANSFER  FUNCTION  OF  THE 
FORWARD  PROPAGATION  OF  DIFFUSE  PHOTON 
DENSITY  WAVES  IN  TURBID  MEDIA  WITH 
AN  EMBEDDED  SPHERICAL  INHOMOGENEITY 


THESIS 

Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 
In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Electrical  Engineering 


Deborah  L.  Lasocki,  B.S.E.E. 
Captain,  USAF 


December  1997 


Approved  for  public  release;  distribution  unlimited 


AFIT/EN /ENG/ 97D-01 


ANALYTIC  TRANSFER  FUNCTION  OF  THE 
FORWARD  PROPAGATION  OF  DIFFUSE  PHOTON 
DENSITY  WAVES  IN  TURBID  MEDIA  WITH 
AN  EMBEDDED  SPHERICAL  INHOMOGENEITY 

Deborah  L.  Lasocki,  B.S.E.E. 

Captain,  USAF 

Approved: 


Dr.  Steven  K.  Rogers  Date 

Committee  Member 


Acknowledgements 

First  and  foremost,  I  would  like  to  thank  God  for  the  wonderfully  complex  world  He  has 
given  us  to  study.  May  we  always  see  the  beauty  He  has  created  around  us. 

I  would  like  to  thank  everyone  who  supported  my  decision  to  take  on  such  a  difficult 
topic,  in  particular,  my  committee  members  Dr.  Byron  Welsh,  Dr.  Mike  Roggemann,  and 
Dr.  Steve  “Captain  Amerika”  Rogers.  You  all  took  the  time  to  help  when  I  needed  it. 
Many  words  of  thanks  need  to  be  bestowed  upon  Major  Pete  Collins,  my  thesis  advisor, 
for  always  believing  in  me,  even  when  the  research  looked  quite  bleak.  It  was  your  allowing 
me  the  freedom  to  work  at  my  own  pace  and  style  that  I  appreciated  most.  Thank  you. 

I  would  also  like  to  thank  my  sponsor,  Dr.  Chuck  Matson  of  the  Advanced  Imaging 
Directorate  at  Phillips  Labs.  Your  excitement  and  brilliant  dabbling  in  this  field  made 
this  thesis  not  only  possible,  but  enjoyable.  Thank  you  for  your  advice,  guidance,  and 
hard-to-answer  questions  that  made  me  think.  But  most  of  all,  thank  you  for  inspiring  the 
budding  theoretician  in  me  to  continue  challenging  the  world  around  me. 

Also  my  parents  must  be  recognized  for  instilling  in  me  the  need  for  higher  education. 
Your  support  helped  me  to  pursue  and  complete  this  thesis.  I  am  only  now  beginning  to 
realize  how  much  your  guidance  throughout  my  life  has  contributed  to  my  current  success. 

A  special  word  of  thanks  to  Christine  Kerschus  for  your  unfailing  optimism,  prayers, 
and  much  needed  laughter  and  study  breaks. 

Finally,  I  want  to  thank  my  future  husband,  Troy  Van  Caster.  You  gave  me  endless 
encouragement  and  patience  during  these  last  16  months,  without  which  I  could  not  have 
kept  the  peace  of  mind  to  complete  this  work.  I  only  hope  I  was  able  to  give  you  the  same. 
I  continue  to  thank  God  for  bringing  us  together  here  at  AFIT;  you  are  my  soul  mate  and 
my  best  friend.  I  am  looking  forward  to  many  future  years  of  exploring  the  world  together. 

“It  is  more  important  to  have  beauty  in  one’s  equations  than  to  have  them  fit  experi¬ 
ment. ..because  the  discrepancy  may  be  due  to  minor  features  which  are  not  properly 
taken  into  account  and  which  will  be  cleared  up  with  further  developments  of  the  the¬ 
ory.. ..It  seems  that  if  one  is  working  from  the  point  of  view  of  getting  beauty  into  one’s 
equations,  and  if  one  has  a  really  sound  instinct,  one  is  on  a  sure  line  of  success  [25]” . 

Deborah  L.  Lasocki 


m 


Table  of  Contents 


Page 

Acknowledgements .  iii 

List  of  Figures .  viii 

List  of  Tables .  x 

List  of  Symbols .  xi 

Abstract .  xiii 

I,  Introduction .  1-1 

II.  Diffuse  Photon  Density  Wave  Propagation  in  Turbid  Media .  2-1 

2.1  Introduction .  2-1 

2.2  Diffusion  Approximation  to  the  Transport  Equation .  2-2 

2.3  Diffuse  Photon  Density  Waves .  2-3 

2.4  Physics  Model  of  Scattered  DPDWs .  2-5 

2.4.1  An  Analytic  Solution  .  2-6 

2.4.2  Semi-Infinite  Media .  2-8 

2.4.3  Resolution  Limitations  of  the  Analytic  Solution  .  .  .  2-9 

2.5  Fourier  Optics  Model  of  Scattered  DPDWs .  2-9 

2.5.1  Wave  Propagation  in  Source-Free  Homogeneous  Media  2-10 

2.5.2  Wave  Propagation  in  Heterogeneous  Media .  2-12 

2.5.3  Reconstruction  Algorithm .  2-13 

2.5.4  Simulated  Results .  2-14 

2.5.5  Experimental  Results .  2-15 


IV 


Page 

III.  An  Analytic  Transfer  Function .  3-1 

3.1  Introduction .  3-1 

3.2  Transfer  Function  Formulation .  3-1 

3.3  Decomposition  of  a  Spherical  Wave  Into  Plane  Waves  ....  3-4 

3.3.1  Incident  Field  Representation .  3-4 

3.3.2  Scattered  Field  Reformulation  .  3-6 

3.4  Fourier  Transform  of  the  Analytic  Solution .  3-9 

3.4.1  Incident  Wave  .  3-9 

3.4.2  Scattered  Wave .  3-10 

3.4.3  Transfer  Function  .  3-10 

3.5  Summary .  3-11 

IV.  Implementation  of  Transfer  Function  Models  and  Test  Plan .  4-1 

4.1  Introduction .  4-1 

4.2  System  Geometry .  4-1 

4.3  Approximations  in  Simulations .  4-2 

4.3.1  Series  Truncation .  4-2 

4.3.2  Sampling  Around  Singularities .  4-3 

4.4  Correlation  Methodology .  4-4 

4.5  Implementation  of  Analytic  Solution  Method .  4-4 

4.6  Analytic  Solution  Simulation  Code  Development .  4-5 

4.6.1  Code  Structure .  4-5 

4.6.2  Input  Parameters .  4-5 

4.6.3  Output  Quantities .  4-6 

4.6.4  Validation  of  Simulation .  4-6 

4.7  Implementation  of  Analytic  Transfer  Function  Method  ....  4-8 

4.8  Analytic  Transfer  Function  Simulation  Code  Development  .  .  4-8 

4.8.1  Code  Structure .  4-9 


v 


Page 

4.8.2  Discretization  Effects  .  4-9 

4.8.3  Input  Parameters .  4-9 

4.8.4  Output  Quantities .  4-9 

4.8.5  Validation  of  Simulation .  4-10 

4.9  Implementation  of  the  Fourier  Optics  Model .  4-11 

4.9.1  Code  Structure .  4-12 

4.9.2  Discretization  Effects  .  4-13 

4.9.3  Input  Parameters .  4-13 

4.9.4  Output  Quantities .  4-14 

4.10  Test  Plan .  4-14 

4.10.1  Perturbation  Cases .  4-14 

4.10.2  Analysis  Plan .  4-15 

4.11  Summary .  4-16 

V.  Simulation  Results  and  Analysis .  5-1 

5.1  Introduction .  5-1 

5.2  Sensitivity  to  Contrasts  in  Optical  Parameters .  5-1 

5.2.1  System  Configuration .  5-2 

5.2.2  Absorption  Contrast  With  a  Pure  Absorber .  5-3 

5.2.3  Scattering  Contrast  With  a  Pure  Scatterer  .  5-5 

5.2.4  Low  Contrast  in  Absorption  and  Varied  Contrast  in 

Scattering  Coefficients .  5-5 

5.2.5  High  Contrast  in  Absorption  and  Varied  Contrast  in 

Scattering  Coefficients .  5-9 

5.2.6  Low  Contrast  in  Scattering  and  Varied  Contrast  in 

Absorption  Coefficients .  5-11 

5.2.7  High  Contrast  in  Scattering  and  Varied  Contrast  in 

Absorption  Coefficients .  5-14 

5.2.8  Summary  of  Sensitivity  Analysis .  5-14 


vi 


Page 

5.3  Sensitivity  to  Inhomogeneity  Size .  5-19 

5.4  Scattering  Moment  Contributions .  5-19 

5.4.1  Pure  Absorption .  5-21 

5.4.2  Pure  Scatterer .  5-22 

5.4.3  High  Absorption  Contrast .  5-22 

5.4.4  Low  Absorptive  Contrast .  5-24 

5.4.5  High  Scattering  Contrast .  5-24 

5.4.6  Low  Scattering  Contrast .  5-25 

5.4.7  Summary  of  Moment  Analysis .  5-26 

5.5  Fourier  Optics  Model  Validity  .  5-35 

5.6  Summary .  5-46 


VI.  Summary .  6-1 

6.1  Conclusions .  6-1 

6.2  Recommendations  for  Future  Research .  6-2 

6.3  Closing .  6-3 

Appendix  A.  Input  Parameter  Code  Listing  to  the  Analytic  Solution  .  .  .  A-l 

Appendix  B.  Analytic  Solution  Code  Listing .  B-l 

Appendix  C.  Transfer  Function  Code  Listing  to  the  Analytic  Solution  .  .  C-l 

Appendix  D.  Fourier  Optics  Wave  Propagation  Code  Listing .  D-l 

Bibliography .  BIB-1 

Vita .  VITA-1 


List  of  Figures 


Figure  Page 

2.1.  Coordinates  for  Inhomogeneity  Geometry .  2-4 

2.2.  Diffuse  Photon  Density  Wave  Propagation  Representation .  2-5 

2.3.  Extrapolated  Surface  for  Semi-Infinite  Medium .  2-9 

2.4.  Fourier  Optics  Reconstruction  Schematic .  2-13 

3.1.  Linear  System  Representation .  3-1 

3.2.  Analytic  and  Fourier  Optics  Solution  Systems .  3-3 

3.3.  Path  of  Integration  in  Transform  Space .  3-7 

3.4.  Region  of  Validity  of  the  Multipole  Expansion .  3-9 

4.1.  System  Geometry  for  Simulations .  4-2 

4.2.  Series  Truncation  Effects  on  Precision .  4-3 

4.3.  Percent  Error  of  Analytic  Incident  Wave .  4-8 

4.4.  Percent  Error  of  Analytic  Scattered  Wave .  4-9 

4.5.  Percent  Error  of  Analytic  Total  Wave .  4-10 

4.6.  Flow  Chart  for  Analytic  Transfer  Function  Method .  4-11 

4.7.  Percent  Error  of  Analytic  Transfer  Total  Wave .  4-12 

4.8.  Flow  Chart  for  the  Fourier  Optics  Model .  4-13 

5.1.  Location  of  Maximum  Perturbation  in  Lateral  Direction .  5-2 

5.2.  Homogeneous,  Scattered,  and  Total  Wave  Slices  of  a  Pure  Absorber  5-4 

5.3.  Analytic  Transfer  Function  for  Pure  Absorbers .  5-6 

5.4.  Homogeneous,  Scattered,  and  Total  Wave  Slices  of  a  Pure  Scatterer  5-7 

5.5.  Analytic  Transfer  Function  for  Pure  Scatterers .  5-8 

5.6.  Analytic  Transfer  Function  for  Low  Contrast  in  /ia  .  5-10 

5.7.  Analytic  Transfer  Function  for  High  Contrast  in  (ia .  5-12 


VUl 


Figure  Page 

5.8.  Analytic  Transfer  Function  for  Very  High  Contrast  in  /xQ .  5-13 

5.9.  Analytic  Transfer  Function  for  Low  Contrast  in  //' .  5-15 

5.10.  Analytic  Transfer  Function  for  High  Contrast  in  n's .  5-16 

5.11.  Analytic  Transfer  Function  for  Very  High  Contrast  in  fi's .  5-17 

5.12.  Sensitivity  to  Change  in  Inhomogeneity  Size .  5-20 

5.13.  Moments  of  a  Pure  Absorber,  High  Contrast .  5-27 

5.14.  Moments  of  a  Pure  Absorber,  Low  Contrast .  5-28 

5.15.  Moments  of  a  Pure  Scatterer,  High  Contrast .  5-29 

5.16.  Moments  of  a  Pure  Scatterer,  Low  Contrast .  5-30 

5.17.  Moments  for  High  Contrast  in  fia  and  fi's .  5-31 

5.18.  Moments  for  High  Contrast  in  fia  and  Low  Contrast  in  fi's  .  5-32 

5.19.  Moments  for  Low  Contrast  in  fia  and  /i'  5-33 

5.20.  Moments  for  Low  Contrast  in  fia  and  High  Contrast  in  fi's  .  5-34 

5.21.  Mean  Square  Error  for  Pure  Absorbers .  5-38 

5.22.  Mean  Square  Error  for  Pure  Scatterers .  5-39 

5.23.  Mean  Square  Error  for  Low  Contrast  in  na .  5-40 

5.24.  Mean  Square  Error  for  High  Contrast  in  fj,a .  5-41 

5.25.  Mean  Square  Error  for  Very  High  Contrast  in  /zc .  5-42 

5.26.  Mean  Square  Error  for  Low  Contrast  in  n's .  5-43 

5.27.  Mean  Square  Error  for  High  Contrast  in  n's .  5-44 

5.28.  Mean  Square  Error  for  Very  High  Contrast  in  /j,'s .  5-45 


IX 


List  of  Tables 

Table  Page 

4.1.  Mean-Square  Error  Verification  of  Simulation  Codes . 4-12 

4.2.  Typical  Optical  Parameters  of  Tissue .  4-14 

4.3.  Optical  Parameter  Test  Matrix .  4-15 


x 


List  of  Symbols 


Symbol  Page 

L  spectral  radiance  .  2-2 

Cl  direction .  2-2 

t  time .  2-2 

v  speed  of  light  in  medium .  2-2 

Ht  transport  coefficient  .  2-2 

fia  absorption  coefficient  .  2-2 

Us  scattering  coefficient .  2-2 

y's  reduced  scattering  coefficient .  2-2 

gi  scattering  angle .  2-2 

f(Cl,  Cl')  normalized  phase  function  .  2-2 

S(f,  Cl,  t)  angular  distribution  of  the  source .  2-2 

D  diffusion  coefficient .  2-3 

So(f,t)  isotropic  moment  of  the  source  .  2-3 

oj  angular  frequency .  2-4 

k  complex  wavenumber .  2-4 

Sac  source  modulation  amplitude .  2-4 

$AC  homogeneous  wave  solution .  2-4 

$  total  wave  solution .  2-6 

<E>SC  scattered  wave  solution .  2-6 

jl  spherical  Bessel  function  of  the  first  kind,  order  l .  2-6 

n\  spherical  Neumann  function  of  the  first  kind,  order  l .  2-6 

Y™  spherical  harmonics  of  degree  l,  order  m .  2-6 

Ai>m  scattered  wave  amplitude .  2-7 

zs  source  distance .  2-7 

a  inhomogeneity  radius .  2-7 


xi 


Symbol  Page 

D out  photon  diffusion  coefficient  of  background  medium .  2-7 

D in  diffusion  coefficient  inside  the  inhomogeneity .  2-7 

kout  wavenumber  of  background  medium .  2-7 

kin  wavenumber  inside  the  inhomogeneity .  2-7 

h!jp  spherical  Hankel  function  of  the  first  kind,  type  l .  2-7 

hp)1  first  derivative  of  the  spherical  Hankel  function  of  the  first  kind,  type  l  2-7 

j[  first  derivative  of  the  spherical  Bessel  function  of  the  first  kind .  2-7 

zj  measurement  (detector)  plane  distance .  2-9 

T  Fourier  transform  operator .  2-11 

T~x  Inverse  Fourier  transform  operator .  2-11 

fx  spatial  frequency  variable .  2-11 

fy  spatial  frequency  variable .  2-11 

t/ 1  transmittance  function .  2-12 

Az  2-plane  interval .  2-13 

£  complex  frequency  variable .  3-4 

r)  complex  frequency  variable .  3-4 

C  complex  frequency  variable .  3-4 

Pi  Legendre  polynomial .  3-6 

a  complex  spherical  angle .  3-6 

/3  real  spherical  angle .  3-6 

Uac  homogeneous  wave  solution  in  Fourier  space  .  3-10 

Use  scattered  wave  solution  in  Fourier  space .  3-10 

cr2  mean-square  error .  4-4 


xn 


AFIT/EN /EN  G  /  97D-01 


Abstract 

Diffusing  photons  can  be  used  to  detect  and  localize  optical  inhomogeneities  em¬ 
bedded  in  turbid  media  such  as  clouds,  fog,  paint  and  human  tissue.  This  thesis  shows 
that  a  transfer  function  derived  from  an  analytic  solution  of  the  Helmholtz  equation  can 
completely  characterize  in  three  dimensions  the  perturbations  in  the  forward  propagation 
phenomena  caused  by  a  spherical  defect  object  in  a  multiple-scattering  medium. 

Two  models  of  the  forward  propagation  behavior  of  diffuse  photon  density  waves 
(DPDW)  in  homogeneous,  infinite,  turbid  media  that  contains  a  spherical  inhomogeneity 
are  examined.  DPDWs  are  generated  from  sinusoidally,  intensity-modulated  sources  of 
light  and  are  highly  damped  waves  with  a  complex  wavenumber.  These  waves  can  be 
described  by  expressions  derived  from  the  diffusion  approximation  to  the  linear  transport 
equation.  The  first  model  is  an  exact  analytic  solution  based  on  a  modal  expansion  in 
spherical  harmonics.  The  second  model  uses  Fourier  optics  theory  for  wave  propagation 
in  a  plane  through  homogeneous  turbid  media  containing  a  spherical  lens.  The  Fourier 
optics  model  is  found  to  be  a  good  approximation  to  the  exact  analytic  solution  when 
the  optical  absorptive  contrast  of  the  inhomogeneity  and  the  surrounding  media  is  weakly 
perturbative,  and  the  detector  is  not  near  the  inhomogeneity. 

Using  linear  systems  theory,  a  transfer  function  from  the  analytic  model  is  derived. 
This  function  improves  the  Fourier  optics  model  by  replacing  the  spherical  lens  approx¬ 
imation  with  an  exact  representation  of  the  system  perturbation  behavior.  The  transfer 
function  is  shown  through  simulation  to  completely  characterize  the  sensitivity  of  the  sys¬ 
tem  to  detect  and  localize  in  three  dimensions  inhomogeneities  of  varying  optical  contrasts 
with  the  surrounding  media. 
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FORWARD  PROPAGATION  OF  DIFFUSE  PHOTON 
DENSITY  WAVES  IN  TURBID  MEDIA  WITH 
AN  EMBEDDED  SPHERICAL  INHOMOGENEITY 

I.  Introduction 

Concern  over  the  limitation,  cost-effectiveness,  and  risks  associated  with  ionizing 
radiation  in  the  imaging  and  diagnosis  of  human  breast  cancer  is  motivating  researchers  to 
seek  other,  safer  methods  of  detecting  breast  lesions  [6] .  Light  radiation  as  an  imaging  tool 
was  apparent  as  early  as  1929  in  the  pioneering  mammography  studies  of  Cutler  [20,21]. 
Light  provides  a  non-invasive  probe  to  image  benign  and  malignant  structures  within 
human  tissue.  In  spite  of  initial  support  for  the  new  transillumination  technique,  the 
procedure  fell  into  disrepute,  probably  due  to  inadequacy  of  the  available  equipment  [36]. 
Later  improvements  to  Cutler’s  work  reported  that  transmission  of  light  through  human 
skin  is  decreased  in  the  shorter  wavelength  spectrum  and  is  increased  in  the  near  infrared 
[2,15,43].  This  particular  approach  of  illuminating,  or  diffusing,  near  infrared  photons 
(NIR)  in  turbid  media,  such  as  human  tissue,  is  continuing  to  emerge  as  a  potentially 
important  imaging  modality  [47,54,58]. 

There  are  many  methods  in  which  light  is  applied  to  image  breast  lesions.  One 
of  these  methods  is  based  on  the  concept  of  light  traveling  unscattered,  or  ballistically, 
through  turbid  media.  The  unscattered  light  is  then  distinguished  from  the  scattered  light 
using  time  or  spatial  gating  techniques  [8,28,37,47,50,55].  For  many  biological  tissues, 
the  absorption  length  for  NIR  light  is  much  longer  than  the  scattering  length  and  is  much 
smaller  than  the  dimensions  of  the  sample.  So  in  thick  tissue,  this  ballistic  method  is 
unsuitable  since  unscattered  light  is  highly  attenuated  and  becomes  indistinguishable  from 
the  scattered  light.  When  the  turbid  medium  is  thick,  other  methods  for  non-invasive 
imaging  can  be  applied.  The  traveling  photons  through  the  turbid  medium  are  scattered 
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multiple  times  before  they  are  absorbed  or  transmitted  through  the  medium.  This  process 
of  traveling  through  turbid  media  can  be  accurately  described  by  a  diffusional  process 
[12,40,45], 

The  general  behavior  of  photons  propagating  in  an  absorptive  and  scattering  medium 
is  described  by  the  linear  transport  equation  [9,16,22,31].  Analytic  solutions  to  this  equa¬ 
tion  are  typically  difficult  to  handle  and  computationally  intensive.  Under  certain  condi¬ 
tions,  the  transport  equation  can  best  be  approximated  by  the  photon  diffusion  equation. 
When  the  intensity  of  the  source  of  light  is  sinusoidally  modulated,  the  photon  diffusion 
equation  can  be  rewritten  as  the  Helmholtz  equation.  The  modulated  light  source,  referred 
to  as  a  diffuse  photon  density  wave  (DPDW),  produces  a  light  energy  density  which  propa¬ 
gates  outwards  from  the  source  through  the  turbid  medium  [12,28,44].  On  the  microscopic 
level,  photons  walk  in  random  paths,  but  macroscopically,  the  photons  combine  to  produce 
a  scalar  wave  of  light  energy  density  that  maintains  a  well-defined  wavelength,  amplitude, 
and  phase  [58]  and  can  thus  be  treated  within  the  framework  of  wave  phenomena.  The 
wave  behavior  of  DPDWs  have  been  well  researched,  particularly  the  perturbations  caused 
by  optical  inhomogeneities  [11,17,28,44,51,54]. 

Carcinomas  compared  to  healthy  tissue  absorb  light  differently  in  narrow  spectra, 
particularly  in  the  NIR  spectrum  [36].  Lesions  are  then  typically  modeled  as  spherical 
inhomogeneities  which  have  different  optical  (absorptive  and  scattering)  properties  than 
the  breast  tissue  medium  surrounding  it.  When  a  DPDW  is  incident  on  the  tissue,  the 
contrast  in  the  optical  properties  of  the  inhomogeneity  and  the  background  tissue  causes 
distortions  in  the  DPDW.  Measurements  of  these  distortions  can  be  used  to  characterize 
and  localize  these  inhomogeneities  [12].  The  slab  geometry  is  extensively  used  as  a  model  of 
a  physically  compressed  breast  as  it  occurs  during  mammography.  However,  the  geometry 
can  be  modified  as  an  infinite  medium  if  the  tumor  is  small,  since  the  boundaries  of  the 
medium  do  not  significantly  affect  the  DPDW. 

Analytic  solutions  have  been  determined  to  exactly  model  this  forward  propagation  of 
the  DPDW  in  a  turbid  medium  that  contains  an  embedded  spherical  inhomogeneity  [11,17, 
28,44,51,54].  Another  solution  approach  of  this  phenomena  uses  Fourier  optics  theory  [40]. 
When  the  optical  properties  are  weakly  perturbative,  meaning  that  the  difference  between 
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the  optical  properties  inside  and  outside  the  inhomogeneity  is  small,  the  spherical  tumor 
can  be  approximated  as  a  thin  lens.  The  Fourier  optics  approach  models  the  source  wave 
as  a  superposition  of  plane  waves  traveling  through  a  piece-wise  homogeneous  medium. 
However,  as  the  contrast  in  optical  properties  increases,  the  model  of  the  inhomogeneity 
becomes  invalid. 

In  this  research,  an  exact  transfer  function  is  desired  to  model  an  infinite,  homoge¬ 
neous,  turbid  medium  with  an  embedded  spherical  inhomogeneity.  This  transfer  function 
completely  characterizes  the  system  and  can  be  used  to  exactly  model  the  forward  propa¬ 
gation  of  the  DPDW  using  Fourier  optics  theory.  This  model  can  be  used  to  improve  the 
Fourier  optics  approach  by  modeling  perturbations  caused  by  the  inhomogeneity  exactly 
rather  than  by  approximating  the  object  as  a  thin  lens.  In  addition,  the  transfer  function 
can  be  analyzed  directly  to  determine  the  sensitivity  of  the  system  to  detect  and  localize 
different  inhomogeneities,  particularly  those  which  weakly  perturb  the  system. 

This  thesis  is  structured  into  six  chapters  covering  the  history  of  optical  imaging,  the 
development  and  the  simulation  of  the  transfer  function,  and  a  comparative  analysis  of  sev¬ 
eral  forward  propagation  models  used  in  research  today.  Specifically,  Chapter  II  reviews 
the  photon  transport  equation,  the  diffusion  model,  and  the  optics  associated  with  the 
scalar  wave  (DPDW)  solution  that  arise  from  the  model.  In  addition,  aspects  of  Fourier 
optics  theory  are  presented  which  relate  to  its  application  to  the  forward  propagation  of  a 
spherical  wave  through  a  heterogeneous  medium.  Chapter  III  develops  an  exact  transfer 
function  to  completely  characterize  a  spherical  wave  traveling  through  an  infinite  homo¬ 
geneous  medium  containing  an  embedded  spherical  inhomogeneity.  Chapter  IV  discusses 
the  structure  of  the  simulations,  demonstrates  the  validity  and  accuracy  of  the  transfer 
function  method  as  well  as  the  Fourier  optics  model,  and  outlines  the  specific  inhomo¬ 
geneity  optical  characteristics  that  are  of  interest  in  this  thesis.  Simulation  codes  for  the 
transfer  function  and  the  Fourier  optics  methods  are  presented  in  Appendices  C  and  D, 
respectively.  Chapter  V  presents  and  discusses  the  results  of  the  simulations  conducted. 
Chapter  VI  reviews  conclusions  made  throughout  the  thesis  and  recommends  areas  for 
further  study. 
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II.  Diffuse  Photon  Density  Wave  Propagation  in  Turbid  Media 
2.1  Introduction 

Visually  opaque  media  are  unusual  in  that  some  are  considered  opaque  because  they 
strongly  absorb  visible  light.  Others  such  as  paint,  foam,  milk,  and  human  tissue  are 
optically  opaque  because  the  light  traveling  within  them  is  highly  scattered.  In  fact,  a 
very  small  number  of  photons  travel  in  straight  lines  through  such  materials.  Instead, 
the  photons  incident  within  a  thick  material  multiply  scatter  and  trace  out  random  paths 
until  they  are  either  absorbed,  or  they  escape  through  the  boundaries  of  the  media.  Media 
that  cause  photons  to  propagate,  or  “diffuse”,  through  the  material  in  that  manner  are 
also  referred  to  as  optically  “turbid”  media.  Human  tissue  acts  as  a  turbid  medium  in  a 
narrow  spectral  window  of  about  700-900  nm  where  the  light  is  mostly  scattered,  and  the 
absorption  effects  are  small. 

The  high  degree  of  scattering,  a  result  of  a  very  short  scattering  length  characteristic 
of  photon  transport  in  turbid  media,  has  been  shown  to  be  very  well  approximated  as 
a  diffusive  process.  This  diffusion  transport  behavior  of  photons  is  commonly  known  as 
photon  migration.  When  the  intensity  of  the  incident  light  wave  is  sinusoidally  modulated, 
it  has  a  small  but  measurable  traveling  wave  disturbance  of  the  light  energy  density  prop¬ 
agating  outwards  from  the  source.  This  is  hence  referred  to  as  a  diffuse  photon  density 
wave  (DPDW)  [12,28,44]. 

In  this  chapter,  photon  migration  properties  are  investigated  through  the  radia¬ 
tive  transport  equation  and  its  approximation  leading  to  the  photon  diffusion  equation. 
DPDWs  are  known  to  adhere  to  standard  scattering  models.  A  particular  analytic  model 
for  localizing  a  spherical  heterogeneity  in  an  otherwise  homogeneous  medium  as  well  as  its 
resolution  limitations  are  presented.  Finally,  a  Fourier  optics  model  to  the  same  perturba¬ 
tion  problem  and  its  resolution  performance  are  discussed.  These  solution  developments 
form  a  basis  for  the  comparison  of  two  approaches  to  deep  tissue  imaging  with  near  infrared 
light. 
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2.2  Diffusion  Approximation  to  the  Transport  Equation 

The  propagation  of  waves  in  scattering  media  can  be  rigorously  described  with  the 
linear  transport  equation  [9,16,22,31] 

-  dL(-r’^  +  v  •  L(r ,  Cl,  t)Cl  +  ptL(f,  n,  t)  = 

v  at 

ps  J  L(r,  Cl',  t)f(Cl,  Cl')dCl'  +  S(f,  Cl,  t ).  (2.1) 

The  gradient  symbol  V  =  ^ i  +  +  Jjfc.  where  i,  j  and  k  are  unit  vectors  in  the  x,  y,  and 

z  directions,  respectively.  L(f,  Cl,  t )  is  the  spectral  radiance,  or  the  number  of  photons  per 
unit  volume,  per  unit  wavelength,  per  unit  solid  angle  (sr  =  steradian  =  unit  solid  angle)  at 
f,  traveling  in  direction  Cl,  at  time  t,  with  units  of  Wm-2sr-1.  v  is  the  speed  of  light  in  the 
transporting  (turbid)  medium  which  is  characterized  by  the  transport  coefficient  m  where 
pt  —  ys+na.  The  absorption  coefficient  is  pa  in  units  of  m-1,  and  the  scattering  coefficient, 
ps,  in  the  same  units,  is  the  reciprocal  of  the  scattering  length.  The  reduced  scattering 
coefficient  p's  is  the  reciprocal  of  the  random  walk  step,  meaning  it  is  the  reciprocal  of  the 
average  path  length  before  the  photon’s  direction  becomes  random.  The  two  scattering 
parameters  are  related  by  the  measure  of  the  amount  of  incident  light  scattered  in  the 
forward  direction,  i.e.  the  average  cosine  of  the  scattering  angle,  g\  [28,35]. 

p's  =  ps(l  -gi)  =  Ms(l-  <  cos0  >)  (2.2) 

The  normalized  phase  function  f(Cl,Cl')  represents  the  probability  of  scattering  from  di¬ 
rection  Cl  into  direction  Cl'.  The  last  term,  S(f,  Cl,  t)  is  the  spatial  and  angular  distribution 
of  the  source  with  units  of  Wm-3sr-1. 

The  transport  equation  is  usually  interpreted  as  a  conservation  equation  for  the  radi¬ 
ance,  L(f,Cl,t).  Many  of  the  difficulties  of  the  solutions  to  the  transport  equation  include 
numerical  analysis  and  heavy  computational  efforts.  The  Pj\r  approximation  method  re¬ 
duces  those  difficulties  and  still  results  in  a  valid  solution  when  applied  to  the  transport 
equation  [10,16,22,31].  The  Pjv  approximation  expands  the  radiance  L,  the  phase  function 
/,  and  the  source  S  into  a  truncated  sum  of  spherical  harmonics  Y^m  [10].  The  Pi  approx- 
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imation  is  valid  if  it  is  assumed  that  the  source-detector  separation  is  large  as  compared 
to  (1  / Us)  ,  and  that  the  scattering  frequency  is  much  larger  than  the  source  modulation 
frequency  [28].  The  condition  fis/i^s  +  Ma)  ~  1  must  also  be  met  for  the  Pi  approxima¬ 
tion  to  be  valid  [30].  When  the  scattering  dominates  the  absorption  in  the  medium,  e.g. 
n's  »  na,  then  the  radiance  can  be  expressed  as  an  isotropic  source  [10,29].  In  cases 
where  photon  absorption  is  large  as  in  the  liver,  hematomas,  and  regions  with  large  blood 
concentrations,  higher  order  solutions  to  the  transport  equation  would  be  a  more  accurate 
approximation.  The  optical  properties  of  the  human  brain  or  breast  demonstrate  that  the 
Pi  approximation  is  valid  [27] .  The  photons  are  treated  as  particles  in  this  work,  thus  the 
interference  terms  from  superposition  effects  as  well  as  polarization  are  neglected  in  this 
simplification  of  the  transport  equation  [26] .  If  the  source  is  sinusoidally  modulated,  then 
applying  these  assumptions  reduces  the  transport  equation  (2.1)  to  the  photon  diffusion 
equation  [10,16] 

-DV2$(f ,  t )  +  vfj,a$(r,  t)  +  ^  =  vS0(f,  t).  (2.3) 

ot 

Here,  the  photon  fluence  (or  energy  density)  is  represented  as 

$(f,i)  =  J  (2.4) 

and  has  units  of  Wm~2.  The  diffusion  coefficient  is  given  by  D=  v/[3([ia  +  n's)]  is  in  units 
of  m-1,  and  So(f,  t)  is  the  monopole  (isotropic)  moment  of  the  source. 

It  has  been  demonstrated  that  at  a  sharp  boundary  for  a  small  random  walk  step 
(l//r' ),  the  transport  equation  better  approximates  the  photon  migration  across  the  bound¬ 
ary  [12].  Otherwise,  the  light  transport  in  the  medium  is  well  approximated  by  the  diffusion 
equation. 

2.3  Diffuse  Photon  Density  Waves 

Many  analytic  solutions  to  the  diffusion  equation  have  been  derived  to  model  plane 
wave  sources  in  semi-infinite  media  [33],  as  well  as  on  slab  [16],  spherical  [12,45,52]  and 
other  media  surfaces  [4],  However,  if  the  intensity  of  the  photon  source  is  sinusoidally 
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Figure  2.1  To  solve  the  frequency  dependent  Helmholtz  equation  for  a  spherical  inho¬ 
mogeneity,  it  is  natural  to  use  spherical  coordinates  with  the  origin  at  the 
center  of  the  object.  The  source  is  situated  along  the  z  axis  (0S  =  7 r),  allow¬ 
ing  azimuthal  symmetry  to  be  exploited.  The  radial  distances  for  the  source, 
detector  and  inhomogeneity  are  indicated  in  the  diagram  [10]. 


modulated  at  angular  frequency  cj,  e.g.  5(f)  =  Si>c(r)  +  Sac{x)^xp(— jut),  then  the 
light  energy  density  within  the  turbid  medium  can  be  written  as  a  combination  of  the 
time-dependent  and  time-independent  parts,  i.e.,  $(f)  =  $vc(x)  +  &Ac(r)exp(— jut)  [58]. 

With  an  oscillating  point  source,  Equation  (2.3)  can  be  rewritten  as  the  Helmholtz 
equation  [28]  at  the  field  point  f  with  spherical  polar  coordinates  (r,  6,  </>)  in  the  frequency 
domain  (denoted  with  subscripts  AC)  as 


V2<&Ac(r)  +  k2<S>Ac{r)  =  -^Sac, 

where  k  is  the  complex  photon  density  wavenumber  (with  j  =  yj—  1), 
k2  _  ~vPa  + 


(2.5) 


(2.6) 


and  Sac  is  the  source  modulation  amplitude.  A  time-dependent  solution  to  Equation  (2.5) 
may  be  expressed  in  the  form  [23]  of  the  real  portion  of  t)  =  3?  [<$> Ac{x)exP{~ jut)]. 

From  this  expression,  the  oscillatory  part  of  the  solution  for  a  homogeneous,  infinite  dense 
random  medium  with  a  modulated  point  source  was  determined  to  be  [28] 


^Ac(^)  = 


V^AC  c?fc|r— r'l 

47rD|f  —  f'|  ; 


(2.7) 
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High 


Figure  2.2  The  DPDW  propagates  from  regions  with  a  high  photon  fluence  rate  (or 
energy  density)  to  those  with  lower  rates  (energy  density).  The  density  of  the 
photon  wave  is  indicated  on  the  amplitude  of  the  solid  curve.  The  directions 
of  the  individual  photons  are  depicted  by  the  arrows  [53]. 

where  this  solution  is  referenced  to  the  center  of  the  spherical  inhomogeneity  which  also 
coincides  with  the  origin  of  the  system.  The  unprimed  variable  (f)  correlates  to  the  location 
of  the  detector  and  the  primed  variable  (f1)  refers  to  the  location  of  the  source  as  seen  in 
Figure  2.1.  It  is  customary  to  drop  the  time-dependent  term  for  this  development.  This 
solution  is  a  highly  damped  wave  which  propagates  with  a  single  phase  velocity  and  is 
referred  to  as  a  DPDW.  This  propagation  is  illustrated  in  Figure  2.2  where  the  density 
of  the  photons  is  denoted  by  the  height  of  the  solid  curve,  and  the  individual  photon 
directions  are  depicted  by  the  arrows  [53] . 

Though  the  complex  wavenumber  indicates  that  the  wave  attenuates  very  rapidly 
within  the  medium,  the  wave  possesses  a  well  defined  wavelength,  amplitude  and  phase  at 
all  points  [58].  The  validity  of  the  DPDW  shows  that  the  diffusional  process  of  sinusoidally 
intensity-modulated  light  incident  on  a  turbid  medium  can  be  treated  within  the  framework 
of  wave  phenomena  [10,28], 

2.4  Physics  Model  of  Scattered  DPDWs 

In  the  presence  of  an  inhomogeneity,  the  DPDW  is  distorted  to  a  certain  extent, 
depending  on  both  the  physical  characteristics  as  well  as  the  optical  properties  of  the 
object.  DPDWs  have  been  extensively  researched  and  shown  to  adhere  to  the  standard 
models  of  wave  behavior  such  as  refraction  [44],  diffraction  [11,28],  scattering  from  localized 
heterogeneities  [10],  interference  [17,51],  and  dispersion  [54],  These  properties  of  DPDWs 
have  been  verified  in  both  biological  models  [28]  and  human  breast  studies  [17].  Many  of 
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the  noted  solutions  above  are  only  valid  if  the  inhomogeneity  is  located  between  the  source 
and  detector.  In  addition,  the  object  must  have  specific  optical  properties  relative  to  its 
surrounding  background.  That  is,  the  object  must  be  highly  absorptive  relative  to  the 
background  medium. 

In  this  section,  a  physics  model  of  the  forward  propagation  of  DPDWs  is  investigated. 
In  particular,  an  analytic  solution  to  the  Helmholtz  equation  (Equation  (2.5))  is  discussed 
for  the  same  piecewise  homogeneous  medium  only  now  with  a  spherical  object  embedded 
in  it.  Unlike  other  solutions  [11, 17, 28, 44, 51, 54]  this  solution  is  for  a  highly  scattering 
inhomogeneity  within  another  strongly  scattering  infinite  medium,  and  later  extended  to 
a  semi-infinite  medium  [12]. 

2-4-1  An  Analytic  Solution.  The  total  field,  <h,  in  a  homogeneous  medium  outside 
an  embedded  spherical  inhomogeneity  is  the  summation  of  the  incident  field  {$>ac)  and 
the  scattered  field  (3>sc)  [35],  meaning 

$  =  $AC  +  ^5C)  (2-8) 

where  <&ac  results  from  the  homogeneous  background  medium,  and  $>sc  results  from  the 
inhomogeneity.  The  Helmholtz  equation  (Equation  (2.5))  was  shown  in  Section  2.3  to 
describe  DPDW  propagation  in  piecewise  homogeneous  medium.  With  that  basis  and  ap¬ 
plication  of  appropriate  boundary  conditions,  the  scattered  field  outside  the  inhomogeneity 
has  the  final  form  of  a  modal  solution  to  Equation  (2.5)  involving  spherical  Bessel  (ji)  and 
Neumann  functions  (n;)  of  the  first  kind,  and  spherical  harmonics  of  degree  l  and  order  m 
(XD  [12],  i.e. 

00  l 

^sc{f)  =  [Ahmji{koutr)+ jAijmni{k0Utr)]YlTn(e,4>),  r  >  a.  (2.9) 

1=0  m=-l 
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Here  the  A^m  coefficients  are 


Ai,m  =  ~j  VSA^k°^1)(k0Utzs)Yl*m(n,  0) 

^  Doutkoutaji{kouta)jl(kina )  ~  D inkinajl{kouta)j i(kjna)  ^ 

^out^out^^l^  ^  {Kuto)jl  {kind)  ~  D inkffldhi^  ^  {k0Uf(i)j^k-in(i) 

where  zs  is  equivalent  to  the  source  distance  along  the  z-axis,  a  is  the  radius  of  the  spherical 
inhomogeneity,  and  D out  (D,;n)  is  the  photon  diffusion  coefficient  outside  (inside)  the  sphere, 
and  similarly  kout  ( kin )  is  the  wavenumber  outside  (inside)  the  sphere,  h j1^  is  the  spherical 
Hankel  function  of  the  first  kind,  type  l,  and  h^'  and  j[  are  the  first  derivatives  of  the 
spherical  hankel  and  bessel  functions  respectively.  Applying  azimuthal  symmetry  about 
z  =  0  since  the  spherical  object  is  at  the  origin  leads  to  A^m  =  0  for  m^O.  The  scattered 
field  simplifies  to  [12] 


®sc{r)  =  Y^Ai  [ ji{koutr )  +  jni(koutr )]  y/0) (0, 4>) 
l=o 

OO 

=  r^a-  (2-11) 
1=0 


This  solution  is  simpler  but  similar  to  the  scalar  version  of  Mie  Theory  for  the  scattering 
of  electromagnetic  waves  from  dielectric  spheres  [34]. 

In  the  numerical  calculation  of  the  scattered  wave  (Equation  (2.11)),  the  number  of 
moments  to  include  of  the  infinite  summation  is  important  to  consider.  To  better  under¬ 
stand  the  dependent  parameters  in  the  scattered  wave,  consider  the  first  three  moments 
of  Equation  (2.11),  namely  the  monopole  (l  =  0),  the  dipole  (l  =  1),  and  the  quadrupole 
(l  =  2)  moments  which  are,  respectively, 


Q  II 
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To  leading  order  in  kouta  and  kina,  assuming  \kouta\  «  1  and  \kina\  «  1,  these  moments 
of  the  scattered  wave  evaluate  to  [13]: 


4c°\r)  = 

$5CX)(f)  = 


®sc\r) 


vSAC  ejkr  eikr' 


47 TO3 


-vA  Ha 


87 T^Dout  r  r' 

vSAC  eikreikT' 
87 r2D0Ut  t  r' 


L  3  j  L 
r ,  11 

D0ut  J 

r  11 
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3 

—3cos6Afj,'s 

vSac  ejkr  ejkr'  f,2  ,  i?*  _  _3_' 

87 T2Dout  r  r'  r'  r'2  _ 

(3  cos2  0  —  l) 

X  VSjOUt  +  3A/x's 


(2.15) 


(2.16) 


(2.17) 


The  difference  in  the  absorption  coefficient  of  the  inhomogeneity,  /xa,m5  and  the  absorption 
coefficient  of  the  background,  fJLa,ouU  are  expressed  as  A fia  =  fia>in  -  fia,out ■  Similarly, 
A/x's  =  pi's  in  ~  fJ-'s  out,  where  this  is  the  difference  between  the  reduced  scattering  coefficients 
of  the  inhomogeneity  {n'Sjin)  and  the  background  (/4,<wt)* 

Note  that  the  first  moment,  to  leading  order,  depends  only  on  the  difference  of  the 
absorptive  coefficients.  The  dipole  and  quadrupole  moments,  however,  depend  only  on 
the  difference  of  the  reduced  scattering  coefficients.  Generally,  the  detectable  moments  of 
the  scattered  DPDW  depend  on  the  optical  contrasts  of  the  object  and  the  surrounding 
medium.  These  dependencies  will  be  further  explored  in  Chapters  IV  and  V. 


2. 4-2  Semi-Infinite  Media.  In  medical  imaging,  non-invasive  methods  are  pre¬ 
ferred  and  inherently  involve  boundaries  between  tissue  and  a  non-scattering  medium  such 
as  water  or  air.  Treating  this  system  as  infinite  as  in  Section  2.4.1  will  certainly  lead  to 
a  certain  degree  of  inaccuracy  since  measurements  are  usually  made  by  placing  the  source 
and  detector  against  the  tissue.  This  procedure  creates  the  boundary  between  scattering 
and  non-scattering  media.  Using  the  method  of  image  sources,  this  system  can  be  trans¬ 
formed  into  an  equivalent  infinite  system  which  contains  the  real  and  image  sources  [46]. 
This  transformation  assumes  a  planar  division  and  a  boundary  condition  of  $  =  0  on  an 
extrapolated  zero  boundary  at  a  distance  zq  =  0.7104/ n's  from  the  actual  boundary,  away 


2-8 


Zd  +  z0 


Figure  2.3  The  source  in  the  semi-infinite  medium  can  be  modeled  with  a  negative  source 
image  centered  about  an  extrapolated  surface  situated  at  a  distance  zq  from 
the  actual  boundary  (zd)  away  from  the  turbid  medium  [19]. 


from  the  diffusive  medium  [12].  A  distance  of  zq  =  2/(3  p's)  has  also  been  used  [16,31,33,47]. 
A  complete  listing  of  appropriate  diffusion-based  values  of  Zq  is  presented  by  Aronson  [3] . 

Figure  2.3  illustrates  the  geometry  of  the  extrapolated  surface.  The  z-plane  boundary 
of  the  turbid  medium  is  denoted  by  Zd,  and  the  extrapolated  boundary  is  at  a  distance  Zd  + 
zq.  An  extrapolated  zero  boundary  is  a  good  approximation  for  the  boundary  conditions  in 
the  semi-infinite  case.  Experimental  comparisons  of  the  boundary  model  with  the  infinite 
model  have  demonstrated  that  the  object  sensitivity  has  not  been  affected  [12]. 


2-4-3  Resolution  Limitations  of  the  Analytic  Solution.  An  extensive  analysis 
of  the  ability  of  DPDWs  to  detect  and  localize  inhomogeneities  was  conducted  by  Boas 
et  al.  Their  findings  indicated  that  tumors  as  small  as  0.3  cm  in  diameter  can  be  un¬ 
ambiguously  detected  when  ^a(object)>  3/xa (background),  or  0.4  cm  when  ^'s (object) > 
1.5/x' (background).  If  the  optical  contrast  of  the  tumor  to  the  background  medium  is  not 
known,  then  the  detectable  tumor  diameter  is  on  the  order  of  1.0  cm  or  larger  [13]. 


2.5  Fourier  Optics  Model  of  Scattered  DPDWs 

An  alternative  approach  to  solving  the  forward  propagation  problem  in  the  frequency 
domain  is  to  apply  the  mathematical  tool,  Fourier  analysis.  This  tool  provides  a  means  of 
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describing  linear  phenomena  such  as  that  occurring  in  the  system  of  interest  here.  When 
applying  Fourier  theory  on  an  imaging  system  where  the  information  being  mathematically 
represented  is  of  a  spatial  nature  (e.g.  a  light  amplitude  or  intensity  distribution  over 
space),  the  theory  is  referred  to  as  Fourier  optics  [32]. 

A  new  approach  to  localize  an  object  in  a  turbid  medium  developed  by  Matson  [40] 
uses  Fourier  optics  theory.  The  approach  localizes  a  spherical  inhomogeneity  in  three 
dimensions  using  a  single  set  of  two-dimensional  measurements.  This  solution  method 
includes  several  new  aspects  that  obviates  the  need  for  movement  in  the  source-detector 
setup  or  the  medium  sample  during  measurements,  unlike  solutions  mentioned  in  Sec¬ 
tion  2.4.  Since  only  a  single  view  is  needed,  position  alignment  inaccuracies  are  avoided, 
and  ease-of-use  for  in  vivo  (internal)  imaging  in  a  clinical  environment  is  maintained.  Mul¬ 
tiple  detectors  are  phased  in  a  manner  such  that  the  DPDW  is  scanned  through  the  sample 
providing  the  lateral  localization  data  of  the  inhomogeneity.  Post-processing  techniques 
are  then  used  to  determine  the  depth  information. 

In  this  section,  the  Fourier  optics  theoretical  and  algorithmic  developments  for  the 
forward  propagation  of  a  DPDW  in  homogeneous  media  without  and  with  an  embedded 
scatterer  are  outlined.  The  reconstruction  (or  back-propagation)  algorithm  is  also  pre¬ 
sented  to  show  the  application  of  the  forward  analysis.  Resolution  values  as  found  through 
simulation  are  presented.  In  addition,  current  research  being  conducted  to  validate  this 
novel  approach  is  discussed. 

2.5.1  Wave  Propagation  in  Source-Free  Homogeneous  Media.  An  alternative 
solution  to  describing  the  light  behavior  in  homogeneous  turbid  media  can  be  found  through 
Fourier  optics  due  to  the  linear  nature  of  the  system.  In  Section  2.2,  it  was  shown  that 
light  propagation  in  a  given  medium  can  be  expressed  through  a  solution  to  the  diffusion 
equation  (Equation  (2.3)).  For  a  point  source,  the  result  is  a  spherical  damped  outgoing 
wave.  If  the  medium  is  a  source-free  vacuum,  there  is  neither  scattering  nor  absorption  of 
the  traveling  light,  and  the  solution  to  the  diffusion  equation  is  a  non-damped  wave.  Light 
propagation  in  a  source-free  system  can  similarly  be  analyzed  in  a  linear,  space-invariant 
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(LSI)  system  [32].  Therefore,  using  the  LSI  framework,  Fourier  theory  can  be  applied  to 
determine  the  behavior  of  the  traveling  light. 

In  this  Fourier  optics  approach,  the  Fourier  transform  for  two  independent  variables 
x  and  y  will  be  represented  by  F{g}  and  is  defined  by 

/oo  poo 

/  g(x,y)e~^2n^xX+^yV^dxdy.  (2.18) 

-oo  J — oo 

Similarly,  the  inverse  Fourier  transform  of  a  function  G(fx,fy),  represented  by  JF_1{G}, 
is  defined  as 

/OO  />0 0 

/  G(fx,fy)e^fxX+Mdfxdfy.  (2.19) 

•00  J  —oo 

Since  the  Fourier  transform  is  a  linear  operator,  and  the  incident  and  scattered  fields  are 
linear  systems,  Fourier  transform  properties  can  be  used. 

Consider  an  incident  DPDW  in  the  (x,  y)  plane,  created  by  an  unspecified  system  of 
sources,  and  traveling  in  an  infinite  homogeneous  turbid  medium  in  the  positive  z  direction. 
Since  the  DPDW  satisfies  the  Helmholtz  equation  (Equation  (2.5)),  let  it  be  represented 
by  U(x,y,z  =  0)  across  the  z  =  0  plane  and  its  two-dimensional  Fourier  transform  be 
given  by 

/oo  poo 

/  U(x,y,0)e-^xX+Mdxdy,  (2.20) 

-00  J  —oo 

where  fx  and  fy  are  spatial  frequency  variables.  At  a  parallel  plane  of  positive  distance  z, 
let  Az(fx,  fy)  be  the  two-dimensional  Fourier  transform  of  the  scattered  field  U(x,  y,  z).  In 
other  words,  the  spatial  Fourier  components  are  a  collection  of  various  plane  waves  propa¬ 
gating  in  different  directions  away  from  a  specified  z  plane.  The  amplitude  of  those  Fourier 
components  across  any  other  plane  can  be  identified  as  the  superposition  of  these  plane 
waves  with  the  appropriate  phase  shifts  incurred  during  propagation.  The  relationship  of 
the  Fourier  components  in  different  parallel  planes  is  shown  to  be  [32] 

Az(fx,fy)  =  Ao(fx,  fy;  0)e^V^-(2-A)2-(2-A)2.  (2.21) 
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Within  this  LSI  system,  the  transfer  function,  Hz(fx,  fy),  that  characterizes  the  wave 
propagation  is  [40] 


HZ(fxtfy  \z)  — 


Az(fxi  fy) 

Ao(fxi  fy ;  0) 

^zy/k2- (2nfx)2- (2*7y)2 . 


(2.22) 


In  a  medium  without  embedded  inhomogeneities,  the  scattered  photon  density  wave  can 
be  constructed  through  the  following  relationship  [40]: 


U(x,y;z )  =  T  1 


f  Az(fx,  fy',  zd)  1 

1  Hz{fx,  fy't  z )  J 


(2.23) 


where  Zd  corresponds  to  the  measurement  (or  detector)  plane,  and  2  is  the  location  of  the 
plane  of  interest.  Essentially,  Equation  (2.23)  means  that  the  incident  wave  can  be  mapped 
from  one  z  plane  to  another,  taking  into  account  scattering  and  propagation  effects. 


2.5.2  Wave  Propagation  in  Heterogeneous  Media.  Considering  an  embedded 
spherical  inhomogeneity  in  an  otherwise  piece-wise  homogeneous  medium,  Equation  (2.23) 
is  valid  only  up  to  the  plane  containing  the  object.  To  model  the  diffraction  caused  by 
the  spherical  object,  an  equivalent  dispersive  disk  of  the  same  diameter  can  replace  the 
object  in  the  plane  containing  the  center  of  the  object  [40].  The  ray  optic  model  in 
conventional  optics  is  of  limited  applicability  since  the  measurements  of  the  perturbed 
wave  are  in  the  near  field.  Typically  the  model  only  works  well  for  inhomogeneities  that 
are  not  highly  absorptive  or  dense,  and  have  a  larger  diffusion  coefficient  relative  to  the 
background  medium  [10].  If  the  inhomogeneity  is  replaced  with  this  ray  optic  model,  then 
Equation  (2.23)  is  valid  only  up  to  the  plane  containing  the  center  of  the  object.  The 
diffraction  from  this  disk  is  similar  to  a  thin  lens  with  a  different  wavenumber  value  (kin) 
than  the  surrounding  medium  ( kout ).  The  transmittance  function  associated  with  a  thin 
lens,  tA,  has  been  shown  to  be  [32] 

tA  =  ejAkR,  (2.24) 
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Detection 

Plane 


Figure  2.4  The  photon  density  wave  is  reconstructed  along  planar  slices  (dashed  lines) 
parallel  to  the  measurement  plane.  Each  slice  is  spaced  a  specified  interval 
A 2  [40]. 


where  Ak  is  the  difference  between  the  wavenumber  inside  and  outside  the  inhomogeneity 
(e.g.  Ak  =  kin  ~  kout),  R  is  the  distance  through  the  inhomogeneity  at  (x,y),  specifically 

R  =  2y/a2  -(x-  x')2  -(y-  y')2,  (2.25) 


where  the  primed  and  unprimed  coordinates  are  defined  as  in  Figure  2.1.  The  homogeneous 
wave  is  transmitted  through  the  plane  containing  the  lens  model.  The  remaining  medium 
to  the  right  of  the  inhomogeneity  is  homogeneous,  so  the  now  perturbed  wave  can  be 
propagated  as  in  Equation  (2.21). 

2.5.3  Reconstruction  Algorithm.  To  determine  the  lateral  location  of  the  inho¬ 
mogeneity,  the  photon  density  wave  that  is  due  to  the  background  medium  is  subtracted 
from  the  detected  photon  density  wave  (Equation  (2.8)).  The  lateral  information  in  the 
detection  plane  (Figure  2.4)  can  be  determined  from  either  the  resulting  amplitude  or 
phase  of  the  scattered  photon  density  wave. 

For  depth  information,  the  photon  density  wave  must  be  reconstructed  by  calculating 
planar  slices  of  the  medium  parallel  to  the  measurement  plane  via  Equation  (2.22)  at 
specified  intervals  of  Az.  Equation  (2.23)  is  valid  only  up  to  the  plane  containing  the 
center  of  the  inhomogeneity.  At  that  point,  there  is  a  singularity,  and  the  reconstructed 
photon  density  wave  peaks.  By  determining  where  the  photon  density  wave  peaks  along 
the  z  axis,  the  depth  of  the  inhomogeneity  can  be  found  [40].  Figure  2.4  illustrates  this 
reconstruction  geometry  [40]. 
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Equation  (2.23)  is  the  basic  algorithm,  but  it  is  modified  to  include  a  stabilizing 
factor  for  two  reasons.  First,  the  exponential  fall-off  of  the  transfer  function  causes  the 
reconstructed  Fourier  spectra  to  become  unstable.  Secondly,  when  over  1.0  cm  or  more 
into  the  volume,  computer-roundoff  error  and  any  noise  introduced  into  the  system  be¬ 
come  significant.  A  pillbox  filter  can  be  applied  to  pass  the  unattenuated  reconstructed 
Fourier  spectra  within  a  limited  radial  distance  around  the  origin  to  attain  system  sta¬ 
bility.  Another  correction  is  made  to  the  algorithm  since  exponential  attenuation  of  the 
DPDW  within  the  turbid  media  limits  the  dynamic  range  of  the  image  intensities  during 
reconstruction.  To  remain  within  the  dynamic  range  of  the  reconstructed  image  the  recon¬ 
structed  wave  intensity  is  normalized  by  the  energy  in  the  wave  at  the  plane  of  interest  [40] . 

2.5.4  Simulated  Results.  Using  the  algorithm  described  in  the  previous  section, 
computer  simulation  results  [40]  demonstrated  the  validity  of  the  methodology  to  localize 
a  spherical  object  embedded  in  an  infinite  homogeneous  medium.  With  mammography  in 
mind,  the  medium  was  modeled  after  human  breast  tissue  with  a  spherical  tumor  located 
at  multiple  positions.  The  surrounding  background  breast  tissue  was  characterized  at  700 
nm  light  with  a  reduced  scattering  and  absorption  coefficients  of  14.0  and  0.035  cm-1, 
respectively  [58].  The  tumor  is  modeled  with  carcinoma  tissue  parameters  characterized 
at  a  wavelength  of  780  nm  and  had  reduced  scattering  and  absorption  coefficients  of  12.0 
and  0.5  cm-1,  respectively.  The  radius  of  the  tumor  was  set  at  0.5  cm  [48].  The  source 
was  modeled  as  a  point  source  with  an  amplitude  modulation  of  1  GHz,  giving  a  DPDW 
wavelength  of  2.67  cm.  The  source’s  closest  point  in  the  detection  plane  was  5  cm.  Zero- 
mean  complex  Gaussian  noise  with  1%  amplitude  and  1°  phase  standard  deviations  was 
added  to  computer-simulated  data  generated  with  a  previously  developed  code  [12,40,49]. 

Under  these  simulated  conditions,  the  lateral  localization  of  the  1.0  cm  diameter 
tumor  was  accurate  but  dependent  on  the  amount  of  noise  in  the  system.  Using  the 
algorithm  outlined  in  Section  2.5.3  with  an  interval  distance  A z  =  0.125  cm,  the  depth  of 
the  tumor  was  identifiable  up  to  3.0  cm  into  the  tissue  from  the  detection  plane.  Resolution 
decreased  as  depth  increased,  but  within  a  few  millimeters,  the  center  of  the  tumor  was 
still  determined. 
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2.5.5  Experimental  Results.  Physical  measurements  were  conducted  to  vali¬ 
date  the  algorithm  outlined  in  Section  2.5.3.  The  setup  consisted  of  a  single  amplitude- 
modulated  light  source,  a  multiple  detector  array,  and  a  series  of  tissue  phantoms  (plastic 
resin  with  Ti02  powder  and  NIR  dye)  containing  various  “tumors”  in  the  forms  of  differ¬ 
ent  colored  beads.  The  colors  provided  various  absorption  properties  that  contrasted  with 
the  surrounding  medium.  The  modulation  has  been  lowered  from  1  GHz  in  the  simulated 
experiments  down  to  20-40  MHz  where  preliminary  measurements  showed  more  conclusive 
results.  The  experimental  results  demonstrated  the  ability  of  the  algorithm  to  detect  and 
localize  various  heterogeneities  with  different  optical  properties  relative  to  the  surrounding 
medium  [38]. 
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III.  An  Analytic  Transfer  Function 


3.1  Introduction 

There  are  several  methods  of  determining  the  behavior  of  light  propagation  in  turbid 
media  as  was  shown  in  Chapter  II.  Two  particular  methods  yielded  solutions  that  were 
investigated  in  some  detail. 

The  analytic  solution  solved  the  Helmholtz  equation  for  both  an  incident  wave  in  a 
homogeneous  turbid  medium  and  a  scattered  wave  off  an  inhomogeneity.  This  solution  is 
in  the  form  of  a  multipole  expansion  of  spherical  harmonics. 

An  alternate  method  applied  the  theory  of  Fourier  analysis  to  represent  the  wave  in 
different  regions  of  the  turbid  medium.  In  one  of  the  regions,  a  transfer  function  can  be 
applied  to  propagate  the  perturbed  wave  through  the  media. 

In  this  chapter,  the  objective  is  to  determine  a  single  transfer  function  that  can  be 
applied  to  the  entire  system  by  using  the  analytic  solution  results. 


3.2  Transfer  Function  Formulation 

In  linear  systems,  for  a  given  input  and  output  pair,  there  is  a  transform  that  can 
be  applied  to  the  input  to  yield  the  given  output.  Figure  3.1  illustrates  this  process.  For 
this  transfer  function  development,  the  system  that  is  under  consideration  has  been  shown 
that  it  can  be  modeled  as  a  linear  system  [32],  The  input  to  the  system  is  a  point  source 
of  intensity  modulated  light,  and  the  output  can  be  modeled  in  two  ways.  The  first  being 
in  the  form  of  an  analytic  solution  as  discussed  in  Section  2.4.  The  other  is  the  Fourier 
optics  representation  as  presented  in  Section  2.5. 


X 


H 


Y  =  XH 


Figure  3.1  In  a  linear  system,  for  a  given  input  X,  a  transfer  function  H  can  be  applied 
to  yield  the  output  Y. 
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The  performance  of  both  solutions  have  been  demonstrated  to  detect  and  localize 
inhomogeneities  within  a  turbid  medium  within  a  few  millimeters  [12,40].  Since  the  model 
of  the  inhomogeneity  is  different  in  the  analytic  solution  and  the  Fourier  optics  solution, 
the  systems  are  not  the  same.  The  Fourier  optics  solution  is  applied  piece-wise  to  three 
regions:  the  homogeneous  medium  to  the  left  of  the  inhomogeneity,  the  heterogeneous 
medium  which  contains  the  inhomogeneity,  and  the  homogeneous  medium  to  the  right 
of  the  inhomogeneity.  The  transfer  function  in  the  Fourier  optics  solution  can  only  be 
applied  to  the  homogeneous  regions.  The  analytic  solution  system  encompasses  the  entire 
heterogeneous  medium  comprised  of  both  the  background  and  inhomogeneity.  Thus  the 
analytic  solution  can  be  used  to  determine  a  transfer  function  for  the  entire  system. 

Since  the  total  field  outside  the  inhomogeneity  is  a  superposition  of  the  incident  and 
scattered  waves  (Equation  (2.8)),  then  via  the  linearity  property  of  Fourier  transforms,  the 
Fourier  transform  of  the  total  wave,  represented  by  T  {$},  is  simply  a  superposition  of  the 
individual  Fourier  transforms  of  the  incident  ($^c)  and  scattered  ($sc)  waves.  Meaning, 

.F{$}  =  F{$ac}  +  -7r{$sc}  •  (3.1) 

Rearranging  terms,  Equation  (3.1)  can  be  rewritten  as 

yW=^{$AC}{l  +  ^||M},  (3.2) 

provided  l/T  {$ac}  exists.  The  right  bracketed  term  in  Equation  (3.2)  is  the  transfer 
function.  The  analytic  form  of  the  transfer  function  is  the  ultimate  objective  of  this  chapter. 
The  Fourier  transforms  for  the  incident  and  the  scattered  waves  will  be  developed  and  the 
transfer  function  then  determined. 

Throughout  these  developments,  the  Fourier  transform  for  two  independent  variables 
x  and  y  will  be  represented  by  J7  {g}  and  is  defined  by 

/oo  poo 

/  g{x,y)e~^x+r)y'>dxdy.  (3.3) 

■OO  J  —oo 
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Analytic  Model 


Hl(q) 


Ul(q)=X1(q) 


Ui{q)  =  Xl{q)Hl{q) 


Fourier  Optics  Model 


ui  (z)  =  x2(z) 


H2(q) 


u2(z)  =  X2  (z)tA 


U2(q)  =  X{x2(z)tA}H2(q) 


Figure  3.2  The  top  system  depicts  the  analytic  solution  model  as  a  linear  system  rep¬ 
resented  by  an  input,  Xi(q),  propagated  through  the  medium  containing  a 
spherical  inhomogeneity  by  a  transform  function,  Hi(q),  to  become  the  out¬ 
put,  Ui(q).  The  Fourier  optics  model  in  the  bottom  system  shows  an  incident 
wave  through  homogeneous  medium  and  a  thin  lens  approximation  of  a  spher¬ 
ical  inhomogeneity  (£4),  and  propagated  out  of  the  medium  by  H2{q).  The 
two  systems  are  not  equivalent  due  to  the  different  models  of  the  spherical 
inhomogeneity. 

Similarly,  the  inverse  Fourier  transform  of  a  function  G(£,  77),  represented  by  X~l  {(?},  is 
defined  as1 


Jr_1  {«>  =  d{d„. 

(27TJ  J —00  J — 00 


1This  transform  definition  differs  slightly  from  that  defined  in  Chapter  II.  The  difference  is  in  a  scale 
factor  which  does  not  change  the  properties  of  the  transform. 
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3.3  Decomposition  of  a  Spherical  Wave  Into  Plane  Waves 

The  difficulty  of  the  problem  of  a  spherical  wave  on  a  plane  interface  is  due  to  the 
difference  between  the  symmetries  of  the  wave  and  the  form  of  the  boundary.  While  the 
incident  wave  caused  by  a  point  source  has  spherical  symmetry,  the  boundary  (detector)  is 
a  plane.  It  is  natural  to  express  the  spherical  wave  in  terms  of  plane  waves  in  the  desired 
Cartesian  coordinate  system  in  order  to  evaluate  the  wave  on  the  planar  boundary. 

Consider  an  optically  dense  homogeneous  medium  with  a  spherical  inhomogeneity 
composed  of  a  highly  scattering  medium  and  whose  center  is  located  at  the  origin  of  a 
Cartesian  coordinate  system  ( x,y,z ).  A  diffuse  photon  density  wave  (DPDW)  is  assumed 
incident  in  the  direction  of  the  +z  axis  and  is  observed  at  a  detector  position  in  the  region 
bounded  by  z  >  a,  where  a  is  the  radius  of  the  inhomogeneity. 

In  this  section,  the  plane  wave  formulation  of  the  incident  wave  and  the  scattered 
wave  from  the  spherical  object  are  investigated. 


3.3.1  Incident  Field  Representation.  To  obtain  the  incident  wave  in  Equa¬ 
tion  (2.7),  the  scalar  Green’s  function  was  used  to  solve  the  Helmholtz  equation  (Equa¬ 
tion  (2.5)).  The  Green’s  function  has  been  shown  to  have  an  equivalent  triple  integral 
representation  [5,56],  e.g. 


ejk\f-r'\ 

|f  -  f'| 


oo  poo  poo  pM(x-xl)+r]{y-y')+ C(z-V)] 


2?r2  J —oo  J —oo  J — 


£2  +  V2  +  C2  -  k2 


-dCdr/dC. 


(3.5) 


The  unprimed  variable  (f)  indicates  the  distance  to  the  detector  and  the  primed  variable 
(f')  is  the  distance  to  the  source  as  seen  in  Figure  2.1.  The  complex  wave  number  of  the 
medium  is  denoted  by  k.  The  path  of  integration  for  each  of  the  three  complex  variables  £, 
T],  and  (  is  the  corresponding  real  axis.  To  keep  the  frequency  variables  £  and  rj  vigorously 
real,  the  third  variable  £  needs  to  be  bounded  by  a  strip  region  —Q(k)  <  9(£)  <  5(fc) 
in  order  to  keep  the  function  analytic,  where  3r(-)  is  the  imaginary  part  of  the  complex 
number.  This  boundary  assumes  that  there  is  some  dissipation  in  the  medium,  however 
small,  i.e.  that  k  has  a  positive  imaginary  part. 
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Let  7  be  defined  as 


7  =  \/£2  +  r]2  -  k2,  if  +r)2  >  k2, 

7=  -jy/k2  -£2  -rj2,  if  Z2  +  V2<k2.  (3.6) 

This  definition  requires  the  real  part  of  the  square  root  to  be  greater  than  zero.  The 
constraints  on  7  are  valid  in  that  the  right  side  of  Equation  (3.5)  satisfies  the  Helmholtz 
wave  equation  and  gives  the  correct  value  of  the  field  at  z  —  z'  [14].  However,  in  a 
turbid  medium  the  real  part  of  k 2  is  always  negative,  and  the  imaginary  part  is  always 
positive  [12].  This  means  that  the  second  constraint  in  Equation  (3.6)  will  never  occur. 
As  a  result,  the  expression  for  7  can  be  simplified  to 

7  =  3?  (v^2  +V2  ~  (\/£2  +V2  ~  k2)  >  (3-7) 

where  3?(-)  denotes  the  real  part  of  a  complex  number  [39].  Integrating  Equation  (3.5) 
with  respect  to  (,  the  double  integral  representation  results  as  [5] 

e^r  r\  -  _L  f°°  f °°  (3,8) 

|r  — f,|  ^J-ooJ-oo'y 

The  elegant  result  of  Equation  (3.8)  is  a  decomposition  of  a  spherical  wave  into  a  super¬ 
position  of  elementary  plane  waves  in  the  x  and  y  directions  as  indicated  by  the  exponent. 
The  waves  exponentially  attenuate  in  the  two  directions  away  from  the  plane  z  =  z’  which 
contains  the  source. 

In  keeping  the  transform  coordinates  arbitrary  in  Equation  (3.8),  the  incident  field 
in  Equation  (2.7)  can  be  rewritten  as 

$Ac{f,f)  =  r  r  d£dV,  z-z'>  0,  (3.9) 

o7T  L)  J J  _o0  7 

where  v  is  the  speed  of  light  in  the  medium,  Sac  is  the  source  modulation  amplitude,  and 
D  is  the  diffusion  coefficient.  The  system  geometry  in  this  research  dictates  that  the  source 
and  detector  will  be  on  opposite  sides  of  the  inhomogeneity,  and  the  incident  wave  travels 
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in  the  positive  z  direction.  Thus  z  >  z'  in  all  cases,  and  so  the  absolute  value  operator  in 
the  exponent  of  Equation  (3.8)  is  replaced  with  its  argument. 


3.3.2  Scattered  Field  Reformulation.  The  scattered  wave  off  of  a  spherical  inho¬ 
mogeneity  was  determined  to  be  a  multipole  expansion  involving  spherical  functions.  To 
simplify  the  notation,  the  scattered  wave  result  in  Equation  (2.11)  is  rewritten  as 

00 

$sc(r)  =  r>a,  (3-10) 

1=0 

where  the  abbreviation  Il[0\  in  spherical  coordinates  f(r,  6, <j>),  is  defined  as 


nS0) 


(r)  = 


2/  +  1 
47T 


/i|^(A:r)P/(cos  8), 


r  >  a. 


(3.11) 


Here,  is  the  spherical  hankel  function  of  the  first  kind,  type  Z,  and  the  zonal  (or  sec¬ 
toral)  spherical  harmonics  (order  m  =  0)  have  been  replaced  with  the  equivalent  Legendre 
polynomials,  Pi,  and  normalization  factor  [1,7].  The  coefficients  Ai  are  defined  in  Equa¬ 
tion  (2.10).  As  before,  k  is  the  wavenumber  of  the  surrounding  medium.  The  superscript 
(0)  will  be  assumed  to  carry  throughout  the  rest  of  the  derivation.  The  objective  is  to 
rewrite  Equation  (3.11)  into  a  plane  wave  expansion  in  Cartesian  coordinates.  To  that 
end,  Equation  (3.11)  can  be  expressed  in  an  integral  format  as  [42] 


Ik 


— ~7V  ' ~T~ ~  f  d/9  f  ejfc'rPi(cosa)  sinada. 
2ttj1  V  4f  Jo  Jb  ' 


(3.12) 


The  contour  B  is  chosen  to  be  from  e  — >■  —  joo  +  e,  0  <  e  <  |,  as  illustrated  in  Figure  3.3. 
The  vector  k  is  of  length  k  and  has  the  complex  spherical  angle  a  and  real  spherical  angle 
/ 3 .  It  should  be  noted  that  jl  is  the  imaginary  number  y/—l  to  the  power  of  l  and  is  not  to 
be  confused  with  ji  which  is  the  spherical  bessel  function  of  order  l.  Using  the  spherical 
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) 


Figure  3.3  The  path  of  integration  denoted  by  B  in  the  a  plane  corresponds  to  the 
positive  half  of  the  real  axis  in  the  A  plane.  [23,42]. 

coordinate  to  rectilinear  transform  variables,  namely 


r=  y/x‘2  +  y2  +  z2, 

x  =  rsin0cos<f), 

0  =  cos_1(f), 

y  =  rsin0sin(j), 

<t>=  tan-1(|), 

z  =  r  cos  0, 

(3.13) 

the  dot  product  in  the  exponent  in  Equation  (3.12)  is  given  by 

k  •  f  =  hr  [cos  0  cos  a  +  sin  a  cos  /?  sin  0  cos  (j>  +  sin  a  sin  /3  sin  0  sin  <j>\ , 

=  k  [z  cos  o;  +  x  sin  a  cos  (3  +  y  sin  a  sin  f3\ .  (3-14) 

Substituting  Equation  (3.14)  into  Equation  (3.12),  the  result  is 

n i  =  /27rd (3  [  eifc(z cos a+x sin a cos P+v sin a sin P)p( (cos  a)  sin adct. 

2tvj1  V  47t  Jq  Jb 

(3.15) 

From  this  spherical  representation  in  transform  space,  a  conformal  transformation  into 
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polar  coordinates  is  defined  as 


A  =  k  sin  a,  dX  =  k  cos  ada; 


7 


=  -sJX2  -k2  =  -jk 


cos  a, 


(3.16) 


where  the  path  of  integration  in  the  A  plane  is  restricted  along  the  positive  half  of  the  real 
axis,  e.g.  0  <  A  <  oo.  This  path  corresponds  to  mapping  onto  the  curve  in  Figure  3.3 
which  depicts  the  a  plane  [23].  Under  the  assumption  that  there  is  some  absorption  in 
the  medium  (9(&)  >  0),  the  integral  will  converge  along  the  path  of  integration  illustrated 
in  Figure  3.3.  Making  the  conformal  transformation  in  Equation  (3.16),  the  expression  in 
Equation  (3.15)  becomes 

n ,  =  -i-r I*”  d/3  [°°  Ie-7V'(xAcos^Asin/J)Pj  AdA.  (3.17) 
2irjl  jk  V  4?r  Jo  J0  7  \  k  ) 

The  result  in  Equation  (3.17)  is  in  the  same  form  as  the  double  integral  of  Equation  (3.8). 
Let  the  frequency  variables  £  —  A  cos  p  and  ij  =  A  sin  /?,  where  A  =  \/£2  +  V2  and  AdAd/3  = 
d£dr].  Note  that  (  and  77  are  real  transform  variables,  so  the  limits  of  integration  are  now 
over  the  entire  real  axis.  Using  these  transform  variables,  Equation  (3.17)  can  be  rewritten 
as 


n. 


1  1 

2irjl  jk 


1 21  +  1  r°°  r00 

v  47t  j _00  y_Qo 


— e_72e J’(€*+w)  pt 
7 


(3.18) 


Using  the  results  in  Equation  (3.18),  the  scattered  wave  (Equation  (2.11))  can  now 
be  expanded  into  a  superposition  of  ordinary  plane  waves  in  Cartesian  coordinates,  namely, 


- 1  (f) 


\z\  >  a 


(3.19) 


It  must  be  noted  here  that  the  multipole  expansion  of  Equation  (2.11)  and  the  plane  wave 
expansion  in  Equation  (3.19)  have  different  regions  of  validity.  The  multipole  expansion  is 
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Figure  3.4  The  multipole  expansion  is  valid  throughout  the  region  r  >  a,  but  the  plane 
wave  expansion  is  valid  on  \z\  >  a  [23]. 

valid  outside  the  inhomogeneity  on  r  >  a.  The  integral  expression  represents  the  scattered 
DPDW  outside  the  strip  \z\  <  a  [23].  Figure  3.4  illustrates  these  two  regions. 


3-4  Fourier  Transform  of  the  Analytic  Solution 

Since  the  Fourier  transform  is  a  linear  operator,  and  the  incident  and  scattered  fields 
are  linear  systems,  Fourier  transform  properties  can  be  used.  In  this  section,  the  two- 
dimensional  Fourier  transform  in  Section  3.2  is  found  using  the  incident  and  scattered 
wave  expressions  developed  in  the  previous  section. 


3.4-1  Incident  Wave.  For  the  two  dimensional  Fourier  transform  development  of 
the  incident  DPDW,  let  C0(£,iy,z,z')  be  defined  as 


c0(t,v;z,z') 


Ie-7  (*-*0 
7 


z  -  z'  >  0, 


so  that  Equation  (3.9)  becomes 


(3.20) 


*Ac(f,  ?')  =  /"  j°°  C0((,  m  z,  z')eMx+^dt dV,  z  -  z' >  0.  (3.21) 

Note  that  it  is  assumed  that  the  incident  wave  source  location  is  at  (0,0)  in  the  x'-y'  plane. 
In  Equation  (3.21)  the  form  for  a  two-dimensional  Fourier  transform  of  the  incident  wave, 
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UaCi  is  readily  apparent  via  Equation  (3.4),  e.g. 

UAc(Z,mz')  =  ^f-C0(Z,V‘,z,z'),  z-z'>  0.  (3.22) 

It  is  interesting  to  note  that  even  though  the  incident  wave  is  singular  at  the  source 
location,  the  non-zero  imaginary  component  of  k  and  the  real  frequencies,  £  and  rj,  keep  7 
from  becoming  zero. 


3. 4-2  Scattered  Wave.  Similarly  for  the  two-dimensional  Fourier  transform  de¬ 
velopment  of  the  scattered  wave,  let  Cz(£,  rj;  z)  be  defined  as  follows: 


Cz(Z,V,z) 


(3.23) 


Using  the  definition  in  Equation  (3.23)  and  the  results  in  Equation  (3.19),  the  scattered 
wave  can  be  expressed  as 

(3-24) 

From  Equation  (3.4),  the  two-dimensional  Fourier  transform  of  the  scattered  wave  in  Equa¬ 
tion  (3.24),  Use,  is 


Usc(Z,V,z)  =  ^T^\j^TrAlCz^r)''Z^  1*1  -  a’ 

1=0  y 


(3.25) 


where  the  constant  terms  have  been  brought  inside  the  Fourier  integral. 


3.4-3  Transfer  Function.  Using  the  results  in  Equations  (3.22)  and  (3.25),  the 
transfer  function  can  be  determined  from  Equation  (3.2).  Thus,  the  analytic  form  of  the 
transfer  function  is 


...  f.  > v  47tD 

=  1  +  E„tw+. 


21  +  1 
47T 


A1P1 


(?) 


■•7  z' 


\z\  >  a,  z  >  z' 


(3.26) 
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This  transfer  function  indicates  that  the  piece-wise  homogeneous  medium  with  an  embed¬ 
ded  spherical  inhomogeneity  is  characterized  by  the  location  of  the  source  along  the  z  axis 
as  well  as  by  the  Legendre  polynomials  in  the  (£,  77)  plane.  The  phase  is  further  affected 
by  the  1/k,  Ai  and  Pi  factors  since  these  are  complex  values. 

3.5  Summary 

In  this  chapter,  an  analytic  transfer  function  was  determined  for  a  system  of  infinite 
homogeneous  turbid  medium  with  a  spherical  inhomogeneity  situated  at  the  origin.  The 
transfer  function  completely  characterizes  this  system  so  that  an  input  DPDW  can  be 
applied  to  produce  a  corresponding  perturbed  output  wave  at  a  detector  plane.  The 
transfer  function  shows  that  some  of  the  three-dimensional  Fourier  components  of  the 
system  may  be  determined  directly  from  knowing  the  two-dimensional  Fourier  components 
in  the  detector  plane.  It  is  then  possible  to  determine  the  three-dimensional  structure  of 
the  inhomogeneity  with  only  the  measurements  in  the  detector  plane. 
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IV.  Implementation  of  Transfer  Function  Models  and  Test  Plan 
4-1  Introduction 

Two,  among  the  many,  exact  methods  of  determining  the  total  output  wave  in  a 
homogeneous  turbid  media  containing  a  spherical  inhomogeneity  are  simulated  in  this 
chapter.  These  methods  include  direct  calculation  using  Boas’s  results  in  Equation  (2.11), 
and  by  using  Fourier  methods  in  Equation  (3.2).  The  two  codes  developed  to  perform  both 
methods  of  wave  generation  are  referred  to  as  the  “analytic  solution”  for  Boas’s  results 
and  the  “analytic  transfer  function”  for  the  Fourier  optics  approach.  The  total  wave 
outputs  of  each  code  are  correlated  to  determine  the  validity  of  the  simulation  methods. 
In  addition  to  these  exact  models,  Goodman’s  propagation  method  is  also  modeled  [32], 
The  plane  wave  propagation  approach  can  be  accurately  applied  to  the  movement  of  a 
wave  through  homogeneous  media.  The  wave  is  represented  as  a  superposition  of  plane 
waves,  and  its  propagation  through  homogeneous  media  can  be  characterized  by  a  transfer 
function.  However,  in  the  presence  of  an  embedded  spherical  inhomogeneity,  this  approach 
becomes  an  approximation  since  the  system  is  no  longer  homogeneous.  By  comparing 
the  exact  analytic  transfer  function  model  to  the  transfer  function  characterization  by 
Goodman  for  inhomogeneous  media,  the  relative  error  can  be  determined  for  the  plane 
wave  approximation. 

In  this  chapter,  the  system  geometry  is  described  and  the  method  of  wave  comparison 
is  introduced.  The  implementation  of  each  model  of  wave  propagation  is  simulated  and 
validated.  In  addition  to  the  simulation  code  descriptions,  this  chapter  outlines  in  a  test 
plan  the  various  optical  parameters  that  will  be  used  to  analyze  the  behavior  of  the  analytic 
transfer  function  from  the  Fourier  optics  method.  Following  the  test  plan,  the  validity  of 
the  wave  propagation  model  developed  by  Goodman  is  determined. 

4-2  System  Geometry 

The  system  geometry  is  consistently  structured  throughout  all  of  the  simulations  and 
analysis  of  the  models  of  wave  propagation.  The  sample  size  is  5  cm  x  5  cm  x  5  cm  and 
is  modeled  as  an  infinite  homogeneous,  turbid  medium  with  a  point  source  located  along 
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Detection 

Plane 


Figure  4.1  A  spherical  inhomogeneity  is  located  at  the  origin  of  the  system.  A  point 
source  is  incident  at  2  cm  to  the  left  (- z  direction)  and  the  detector  plane 
is  3  cm  to  the  right  ( +z  direction)  of  the  inhomogeneity.  The  source  loca¬ 
tion  nor  the  origin  are  included  in  the  system  analysis  due  to  the  associated 
singularities.  Increment  division  along  all  axes  is  0.125  cm. 

the  —5-axis.  A  spherical  inhomogeneity  is  located  at  the  origin  of  the  system,  2  cm  from 
the  source  and  3  cm  from  the  detector  plane  along  the  2-axis.  The  sampling  interval  along 
both  the  x  and  y- axes  is  0.125  cm.  Along  the  2-axis,  the  sampling  interval  is  also  0.125 
cm,  but  neither  the  origin  of  the  inhomogeneity  or  the  source  location  is  included  in  order 
to  avoid  singularities  in  the  calculation  of  the  scattered  and  incident  waves  respectively. 
This  geometry  is  illustrated  in  Figure  4.1. 


4.3  Approximations  in  Simulations 

4-3.1  Series  Truncation.  To  implement  the  exact  analytic  solution  numerically, 
the  infinite  series  in  Equation  (2.11)  is  truncated  at  a  user-defined  limit.  If  the  conditions 
\kouta\  «  1  and  \kina\  «  1  are  valid,  where  kout  ( kin )  is  the  wavenumber  outside 
(inside)  the  inhomogeneity  with  radius  a,  only  the  first  three  moments  of  the  series  are 
detectable  in  experiments  due  to  signal-to-noise  limitations  [13].  These  three  moments  are 
dependent  on  the  absorptive  and  scattering  properties  of  the  inhomogeneity.  To  leading 
order,  the  first  moment  (monopole)  depends  on  the  difference  of  the  absorption  coefficients 
inside  and  outside  the  inhomogeneity.  The  second  moment  (dipole)  and  the  third  moment 
(quadrupole)  depend  on  the  difference  of  the  reduced  scattering  coefficients  inside  and 
outside  the  inhomogeneity  (see  Section  2.5.2  for  development). 
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Figure  4.2  The  magnitude  of  the  incident  wave  due  to  a  point  source  is  plotted  at  the 
top  of  the  graph.  The  magnitude  of  the  first  six  moments  (/)  in  the  scattered 
wave  are  plotted  below  the  incident  wave  which  gives  better  than  10-6  value 
of  precision  in  the  total  wave.  The  optical  parameters  used  here  are  na,out  — 
0.05  cm4,  fi'SiOUt  =  10  cm-1,  fia,in  =  0.15  cm-1,  p'sin  =  15  cm-1. 

To  determine  the  number  of  terms,  or  moments,  needed,  the  first  six  terms  of  the 
scattered  wave  are  plotted  for  a  1  cm  diameter  spherical  inhomogeneity  with  an  absorption 
coefficient,  pa,  of  0.15  cm-1,  and  a  reduced  scattering  coefficient,  p!s,  of  15  cm-1  surrounded 
by  a  background  tissue  with  fia  =  0.10  cm-1  and  p's  =  10  cm-1.  For  this  example, 
|fcina|=1.2995  and  |&outaj=0.6413.  These  values  indicate  that  more  than  three  moments 
are  needed  in  the  series  for  acceptable  accuracy.  Figure  4.2  shows  that  to  gain  a  precision  of 
10— 6 ,  at  least  the  first  six  terms  in  the  infinite  series  in  Equation  (2.11)  must  be  included. 

4-3.2  Sampling  Around  Singularities.  Since  the  homogeneous  wave  is  undefined 
at  the  location  of  the  point  source,  the  sampling  in  the  spatial  domain  does  not  contain  the 
point  source  location.  The  scattered  wave  described  by  Equation  (2.11)  is  only  valid  outside 
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or  on  the  boundary  of  the  inhomogeneity,  so  the  coordinates  inside  the  inhomogeneity  are 
not  included  in  the  sampling  set.  Sampling  arbitrarily  close  to  those  locations  is  valid 
and  is  useful  in  determining  the  behaviors  of  the  waves  without  causing  numerical  errors 
dealing  with  the  singularities.  Note  that  Boas’s  exact  code  is  valid  for  both  inside  and 
outside  the  inhomogeneity,  but  it  has  the  singularity  at  the  point  source  location  [49]. 

4-4  Correlation  Methodology 

The  mean-square  error  methodology  is  used  to  correlate  waves  calculated  by  different 
codes.  The  error  is  calculated  by  first  finding  the  difference  in  the  magnitude  (or  phase) 
of  each  respective  sampled  pixel  in  the  wave  calculated  by  the  theoretical  code  and  the 
corresponding  pixel  in  the  wave  calculated  by  the  approximated  code.  The  difference 
is  squared  and  divided  by  the  total  number  of  pixels.  The  error,  cr2,  is  mathematically 
expressed  as, 

N  M  /Tjt 

"2  =  EE— 

i=l  1 

where  Ul  is  theoretical  value  and  Ue  is  the  estimated  value  at  the  ith  and  jth  pixel,  and 
N  X  M  is  the  number  of  pixels  in  the  wave  matrix.  For  the  analysis  in  this  thesis,  this 
correlation  equation  is  used  in  the  x-y  plane  (i.e.,  the  detection  plane)  which  is  square.  So 
the  total  number  of  pixels  is  N2. 

The  theoretical  value  will  be  the  output  of  Boas’s  exact  code  [49].  Since  only  the 
first  six  moments  were  used  in  the  MATLAB©  simulations,  the  values  are  accurate  on  the 
order  of  10~9  for  the  magnitude  and  10-6  for  the  phase.  Mean  square  errors  are  expected 
to  be  at  best  on  the  order  of  10~18  and  10-12  for  magnitude  and  phase,  respectively. 

4-5  Implementation  of  Analytic  Solution  Method 

Photon  Migration  Imaging  (PMI)  code  for  this  exact  solution  has  already  been  devel¬ 
oped  by  Boas  and  extensively  analyzed  for  accuracy  [49].  This  exact  code  only  yields  the 
incident  and  total  wave  as  outputs.  The  scattered  wave  can  be  calculated  by  subtracting 
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the  incident  wave  from  the  total  wave  as  per  Equation  (2.8).  However,  at  the  time  the 
analysis  in  this  thesis  was  conducted,  the  moments  that  comprise  the  scattered  wave  could 
not  be  recovered.  The  code  developed  in  this  thesis,  referred  to  as  the  analytic  solution 
code,  can  decompose  the  scattered  wave  into  its  moments.  These  moments  can  then  be 
analyzed  individually.  The  total  perturbed  wave  is  calculated  in  this  code  and  is  correlated 
to  Boas’  code  results  to  validate  the  implementation  of  the  algorithm. 

4-6  Analytic  Solution  Simulation  Code  Development 

The  analytic  solution  in  Equation  (2.11)  is  implemented  in  MATLAB  5.0©  to  take 
advantage  of  the  multi-dimensional  matrix  capabilities  of  the  language.  The  source  code 
for  the  algorithm  is  listed  in  Appendix  B.  This  algorithm  implementation  is  verified  against 
the  PMI  code  [49]. 

4-6.1  Code  Structure.  The  code  is  structured  as  a  function  so  that  the  input 
parameters  can  be  passed  into  the  algorithm  and  only  the  output  quantities  are  returned. 
The  incident  and  scattered  waves  are  calculated  from  Equation  (2.7)  and  the  first  six 
moments  of  Equation  (2.11).  The  total  perturbed  wave  is  simply  the  superposition  of 
those  two  waves  via  Equation  (2.8). 

4-6.2  Input  Parameters.  The  user  defines  the  following  optical  parameters  of  the 
system: 

•  Muain  -  Absorption  coefficient  inside  the  inhomogeneity  (cm-1). 

•  Musin  -  Reduced  scattering  coefficient  inside  the  inhomogeneity  (cm-1). 

•  Muaout  -  Absorption  coefficient  of  the  background  medium  (cm-1). 

•  Musout  -  Reduced  scattering  coefficient  of  the  background  medium  (cm-1). 

•  f  -  Modulation  frequency  (Hz). 

•  idim  -  Length  of  x-  and  y-axes  (cm). 

•  inc  -  Increments  (cm)  along  x-  and  y-axes. 

•  zinc  -  Increments  (cm)  along  z-axis. 


4-5 


•  zdet  -  Location  of  detector  plane  (cm). 

•  Nin  -  Index  of  refraction  of  the  inhomogeneity. 

•  Nout  -  Index  of  refraction  of  the  background  medium. 

The  following  parameter  variables  are  passed  into  the  analytic  solution  algorithm: 

•  numdeg  -  Truncate  series  after  this  number  of  terms. 

•  Sac  -  Source  modulation  amplitude. 

•  a  -  Radius  of  the  inhomogeneity  (cm). 

•  rsx,  rsy,  rsz  -  Location  of  the  source  along  x-,  y-,  and  z-axes  (cm). 

•  x,  y,  z  -  Vectors  containing  sampling  points  along  x-,  y-,  and  z-axes  (cm). 

•  Dout  -  Diffusion  coefficient  of  the  background  medium. 

•  Vout  -  Speed  of  light  in  the  background  medium. 

•  Kout  -  Wavenumber  in  the  background  medium. 

•  Din  -  Diffusion  coefficient  of  the  inhomogeneity. 

•  Kin  -  Wavenumber  of  the  inhomogeneity. 

The  variables  in  bold  are  defined  by  the  user,  while  the  plain-text  ones  are  calculated 
from  the  user-defined  values.  An  example  of  this  listing  is  in  Appendix  A. 

4-6.3  Output  Quantities.  The  outputs  of  the  program  are  the  incident  homoge¬ 
neous  wave  and  the  scattered  wave  due  to  the  inhomogeneity  at  each  incremental  point 
along  all  axes.  The  total  wave  is  also  an  output  quantity  and  is  simply  the  superposition 
of  the  first  two  quantities.  In  the  code,  these  are  identified  by  the  variables:  Uinc,  Us,  and 
Uout,  respectively. 

4.6.4  Validation  of  Simulation.  To  verify  the  validity  of  the  code  developed 
here,  each  wave  of  this  code  at  the  detection  plane  is  correlated  to  the  corresponding 
wave  in  the  exact  solution  code  [49]  by  calculating  the  mean-square  error  between  the 
two  via  Equation  (4.1).  The  system  is  modeled  as  illustrated  in  Figure  4.1  at  20  MHz 
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with  a  spherical  inhomogeneity  with  a  1  cm  diameter.  The  background  medium  has  ya  = 
0.05  cm-1,  fi's  =  10  cm-1,  while  the  inhomogeneity  has  the  properties  f. ia  =  0.15  cm-1, 
y!s  =  15  cm-1. 

The  PMI  code  uses  an  approximated  diffusion  coefficient,  D  =  v/(3y's)  while  the 
analytic  solution  code  used  the  exact  definition  of  the  diffusion  coefficient,  e.g.  D  = 
v/[3(ya  +  fig)].  In  addition,  the  PMI  code  does  not  include  the  v/D out  factor  in  the  Ai)fn 
coefficients  in  Equation  (2.10).  To  correlate  the  two  codes,  these  corrections  were  made  to 
the  analytic  solution  code. 

With  the  amendments  to  the  analytic  solution  code  to  match  the  algorithm  im¬ 
plemented  in  the  PMI  code,  the  magnitude  and  the  phase  of  the  incident  waves  have 
mean-square  errors,  a2  =  7.8695  -10-17  and  a2  =  6.6781  -10-8,  respectively.  The  scattered 
and  total  waves  are  correlated  in  the  same  manner.  The  magnitude  and  the  phase  of  the 
scattered  waves  have  mean-square  errors,  a2  =  2.9876- 10-9  and  a2  =  6.5293- 10-8,  respec¬ 
tively.  A  relative  percentage  error  for  each  respective  point  of  the  incident  and  scattered 
waves  in  the  detector  plane  are  depicted  in  Figures  4.3  and  4.4  to  show  that  the  resulting 
error  is  evenly  distributed  over  the  plane  indicating  that  differences  are  due  to  computer 
numerical  round-off.  However,  the  majority  of  the  error  in  Figure  4.4(a)  is  due  to  the 
truncation  of  the  infinite  series  after  six  terms.  The  PMI  code  truncates  the  series  after 
more  than  20  terms,  so  the  correlation  error  is  consistent  with  this  reason. 

The  total  output  wave  of  the  analytic  solution  code  is  only  used  to  validate  this 
algorithm  against  Boas’s  results.  With  the  exceptions  noted  before,  the  algorithm  is  a 
good  match.  The  magnitude  and  phase  of  the  total  output  waves  have  mean-square  errors 
of  a2  =  2.6027- 10-9  and  a2  =  1.0473- 10-6.  Similar  to  the  incident  and  scattered  waves,  the 
error  is  caused  by  computer  round-off  and  series  truncation.  The  relative  percentage  error 
for  each  respective  x  and  y  coordinate  in  the  detector  plane  is  illustrated  in  Figure  4.5. 
The  numerical  round-off  error  is  proportional  to  1/f,  so  the  error  is  greater  nearer  the 
center  of  plane  and  decreases  as  the  (x,  y)  values  in  the  plane  increase.  The  mean-square 
errors  are  summarized  in  Table  4.1.  These  small  mean-square  error  values  give  validity  to 
the  incident  and  scattered  wave  calculations,  as  well  as  the  total  wave,  at  the  detection 
plane. 
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Figure  4.3  The  percent  error  in  the  detector  plane  is  plotted  for  the  (a)  magnitude 
and  (b)  phase  of  the  incident  waves  caused  by  a  point  source  incident  wave. 
For  both  figures,  na,out  =  0.05  cm-1,  /i's  out  =  10  cm-1,  pa^n  =  0.15  cm-1, 
Pa, in  =  15  cm_1- 


4-7  Implementation  of  Analytic  Transfer  Function  Method 

The  total  wave  can  also  be  calculated  in  the  spatial  frequency  domain  as  determined 
in  Section  3.2.  This  code  calculates  the  incident  homogeneous  wave  as  well  as  the  scattered 
wave  due  to  the  inhomogeneity  exactly  as  executed  in  the  analytic  solution  code.  With 
these  two  wave  quantities,  the  transfer  function  can  be  calculated  and  used  to  generate 
the  total  output  wave.  In  essence,  this  code  validates  the  existence  of  a  transfer  function 
through  the  correlation  of  the  total  wave  generated  in  both  the  spatial  and  the  spatial 
frequency  domains.  In  Chapter  V,  the  output  wave  in  this  analytic  transfer  function  code 
is  used  to  analyze  the  Fourier  optics  wave  propagation  model,  and  the  behavior  due  to 
various  optical  parameters  is  investigated  . 


4-8  Analytic  Transfer  Function  Simulation  Code  Development 

The  analytic  transfer  function  solution  in  Equation  (3.2)  is  also  implemented  in 
MATLAB  5.0©  for  its  multi-dimensional  matrix  capabilities.  The  source  code  for  the 
algorithm  is  listed  in  Appendix  C.  This  algorithm  implementation  is  verified  against  PMI, 
a  known  working  code  [49]. 
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(a)  (b) 


Figure  4.4  The  percent  error  in  the  detector  plane  is  plotted  for  the  (a)  magnitude  and 
(b)  phase  of  the  scattered  waves  from  a  spherical  inhomogeneity  with  a  1  cm 
diameter.  For  both  figures,  pa,out  =  0.05  cm-1,  =  10  cm-1,  p.a^n  = 

0.15  cm-1,  n's  in  =  15  cm-1.  The  error  is  mainly  due  to  the  truncation  of  the 
infinite  series. 

4-8.1  Code  Structure.  This  analytic  transfer  function  code  is  structured  so  that 
input  parameters  can  be  passed  into  the  algorithm  and  only  the  output  quantities  are 
returned.  The  algorithm  determines  the  incident  homogeneous  and  the  scattered  waves  in 
the  spatial  domain  and  then  converts  the  two  quantities  to  the  spatial  frequency  domain. 
The  transfer  function  is  determined  via  Equation  (3.2)  and  is  used  to  construct  the  total 
perturbed  wave.  The  total  wave  is  then  transformed  back  into  the  spatial  domain  where 
it  is  then  correlated  with  the  total  wave  generated  directly  in  the  spatial  domain  from  the 
PMI  code  [49] .  The  implementation  of  the  analytic  transfer  function  solution  is  illustrated 
in  Figure  4.6. 

4-8.2  Discretization  Effects.  Zero-padding  of  the  wave  matrices  avoids  wrap¬ 
around  error  in  the  Fourier  transform  calculations.  The  zero-pad  matrix  size  is  selected  to 
be  a  square  matrix  of  a  power  of  2  in  order  to  take  advantage  of  the  FFT  algorithms. 

4-8.3  Input  Parameters.  The  user-defined  parameters  and  input  variables  are 
the  same  as  listed  in  Section  4.6.2.  An  example  of  this  listing  is  in  Appendix  A. 

4.8.4  Output  Quantities.  The  outputs  of  the  program  are  the  incident  homoge¬ 
neous  wave  and  the  scattered  wave  due  to  the  inhomogeneity  at  each  incremental  point 
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(a)  (b) 


Figure  4.5  The  percent  error  in  the  detector  plane  is  plotted  for  the  (a)  magnitude  and 
(b)  phase  of  the  total  perturbed  wave  from  a  spherical  inhomogeneity  with 
a  1  cm  diameter.  For  both  figures,  fia,out  =  0-05  cm-1,  Ms, out  =  10  cm-1, 
Ma,in  =  0.15  cm-1,  n's  in  =  15  cm-1. 

along  all  axes.  The  transfer  function  that  is  calculated  from  those  two  waves  and  the 
resulting  total  perturbed  wave  are  also  outputs  to  this  code.  All  of  these  quantities  are 
identified  by  the  variables:  Uinc,  Uscatt,  H,  and  UoutH,  respectively. 

4-8.5  Validation  of  Simulation.  To  verify  the  validity  of  the  code  developed,  the 
total  perturbed  wave  calculated  in  the  detection  plane  must  be  transformed  into  the  spatial 
domain,  then  correlated  to  the  wave  calculated  in  the  exact  solution  code  by  Boas  [49]. 
The  correction  factors  discussed  in  Section  4.6.4  are  applied  in  order  to  correlate  to  the 
same  algorithm.  Namely,  the  approximated  diffusion  coefficient,  D  =  v/(3fi's),  is  used, 
and  the  v/D out  factor  in  the  A^m  coefficients  is  removed  from  the  code  to  match  the  PMI 
algorithm.  The  analysis  in  Chapter  V  using  the  transfer  function  calculated  here  does  not 
include  the  v/D out  factor,  but  does  use  the  exact  diffusion  coefficient,  D  =  v/(3[/ia  +  n's}). 

The  perturbed  wave  is  calculated  from  a  1  cm  diameter  spherical  inhomogeneity, 
situated  in  a  homogeneous  medium  as  in  Figure  4.1  at  20  MHz.  The  background  medium 
has  Ma  =  0.05  cm-1,  n's  =  10  cm-1,  while  the  inhomogeneity  has  the  properties  fia  = 
0.15  cm-1,  Ms  =  15  cm-1  (the  same  as  the  previous  validation). 

The  magnitude  and  phase  of  the  total  output  waves  have  mean-square  errors  of 
cr 2  =  2.6027  •  10-9  and  a2  =  1.0473  •  10-6.  The  mean-square  errors  match  the  direct 
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Figure  4.6  The  flow  chart  above  shows  the  code  structure  for  simulating  the  analytic 
transfer  function  method  in  MATLAB  5.0©. 


evaluation  of  the  total  wave  in  the  analytic  solution  code  (Section  4.6.4),  indicating  that 
the  Fourier  transform  algorithm  of  MATLAB  5.0©  does  not  introduce  additional  error 
into  the  program.  The  relative  percentage  error  for  each  respective  x  and  y  coordinate 
in  the  detector  plane  is  illustrated  in  Figure  4.7.  The  incident  and  scattered  waves  were 
correlated  in  Section  4.6.4.  The  mean-square  errors  are  summarized  in  Table  4.1.  These 
high  correlation  values  give  validity  to  the  existence  of  the  transfer  function  since  there  is 
a  good  match  to  the  total  perturbed  wave. 


4-9  Implementation  of  the  Fourier  Optics  Model 

In  Section  2.5,  a  solution  was  discussed  for  wave  propagation  in  a  source-free  ho¬ 
mogeneous  medium  using  Fourier  optics.  Since  the  system  under  investigation  contains  a 
spherical  inhomogeneity  in  a  otherwise  piece- wise  homogeneous  medium,  the  solution  can 
only  be  applied  in  the  regions  that  do  not  contain  the  inhomogeneity.  This  code  calculates 
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Figure  4.7  The  percent  error  in  the  detector  plane  is  plotted  for  the  (a)  magnitude  and 
(b)  phase  of  the  total  perturbed  wave  from  the  analytic  transfer  method, 
from  a  1  cm  diameter  spherical  inhomogeneity.  For  both  figures,  iia,out  = 
0.05  cm-1,  n's  mt  =  10  cm'1,  A*0,*n  =  0.15  cm-1,  fi'sin  =  15  cm"1. 


Wave 

Analytic 

Magnitude 

Analytic 

Phase 

Analytic 

Transfer 

Magnitude 

Analytic 

Transfer 

Phase 

Homogeneous 

7.8695-10-17 

6.6781-10-8 

7.8695-10-17 

6.6781-10-8 

Scattered 

2.9876-10-9 

6.5293-10"8 

2.9876-10-9 

6.5293-10-8 

Total 

2.6027-10-9 

1.0473-10-6 

2.6027-10-9 

1.0473-10-6 

Table  4.1  The  mean-square  error  of  the  simulated  waves  calculated  using  both  the  ana¬ 
lytic  and  the  analytic  transfer  methods.  Both  methods  are  compared  to  the 
exact  solution  developed  by  Boas  [49].  Note  that  the  FFT  algorithm  used 
in  the  analytic  transfer  method  does  not  introduce  additional  error  into  the 
calculation. 

an  incident  wave  and  propagates  it  through  the  homogeneous  medium  up  to  the  plane 
containing  center  of  the  inhomogeneity.  At  that  point,  the  spherical  object  is  modeled  as 
a  thin  lens  with  a  corresponding  transmittance  function  in  the  plane  containing  its  center 
(Section  2.5.2).  The  wave  is  then  transmitted  through  that  plane  by  multiplying  the  inci¬ 
dent  wave  by  the  transmittance  function  in  Equation  (2.24).  This  distorted  wave  is  then 
propagated  to  the  detector  plane  using  Goodman’s  method  in  Equation  (2.22)  [41].  The 
regions  of  validity  for  this  model  are  investigated  in  Chapter  V. 

4-9.1  Code  Structure.  Similar  to  the  previous  two  model  implementations,  this 
method  is  written  in  MATLAB  5.0©.  The  function  structure  enables  the  user  to  pass 
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Load  Input  Parameters 


Transform  Perturbed  wave  into 
spatial  frequency  domain 


Figure  4.8  The  flow  chart  above  shows  the  code  structure  for  simulating  the  Fourier 
optics  model  in  MATLAB  5.0©  [41]. 


parameters  into  the  algorithm  and  returns  only  the  desired  output  quantities.  The  source 
code  for  the  algorithm  is  listed  in  Appendix  D,  and  a  flow  chart  illustrates  the  implemen¬ 
tation  of  the  Fourier  optics  method  in  Figure  4.8. 


4-9.2  Discretization  Effects.  Zero-padding  of  the  wave  matrices  avoids  wrap¬ 
around  error  in  the  Fourier  transform  calculations.  The  zero-pad  matrix  size  is  selected  to 
be  a  square  matrix  of  a  power  of  2  in  order  to  take  advantage  of  the  FFT  algorithms. 


4-9.3  Input  Parameters.  The  user-defined  parameters  and  input  variables  are  the 
same  as  listed  in  Section  4.6.2.  The  variables  inc,  idim,  and  zinc  are  passed  into  the  program 
in  addition  to  those  listed  in  Section  4.6.2.  The  diffusion  coefficient  for  the  inhomogeneity, 
Din,  is  not  needed  in  the  program.  An  example  of  this  listing  is  in  Appendix  A. 


4-9-4  Output  Quantities.  The  output  identified  as  the  variable  “homo”  is  the 
incident  homogeneous  wave  through  the  homogeneous  medium  without  any  inhomogeneity 
present.  The  only  other  output  quantity  is  the  wave  perturbed  by  the  inhomogeneity  and 
is  specified  by  “Udet”. 

4-10  Test  Plan 

Only  the  first  six  moments  will  be  included  in  these  simulations  for  the  reasons 
stated  in  Section  4.3.  The  pure  absorption  and  pure  scattering  cases  have  been  highly 
investigated  [13].  The  concentration  of  this  thesis  is  for  inhomogeneities  that  have  small 
absorptive  and  highly  scattering  properties  as  well  as  low  contrast  as  compared  to  the 
background  medium  properties.  For  completeness,  high  contrast  cases  (i.e.  pure  absorption 
and  pure  scattering)  are  investigated  to  determine  trends. 


Type 

pa  (cm-1) 

/4( cm-1) 

Human  Breast  Tissue 

0.05 

10.0 

Adipose  Tumor 

0.7 

9.0 

Fibroadenoma  Tumor 

0.5 

7.0 

Carcinoma  Tumor 

0.5 

12.0 

Table  4.2  These  are  typical  optical  parameters  of  human  tissue  characterized  at  an  op¬ 
tical  wavelength  of  700  nm  [12,18]. 

4-10.1  Perturbation  Cases.  Typical  pa  and  p!s  parameters  used  to  model  hu¬ 
man  tissue  and  tumors  are  listed  in  Table  4.2.  For  this  analysis,  the  background  medium 
was  consistently  modeled  as  human  tissue  with  pa  =  0.05  cm-1  and  p's  =  10  cm-1.  The 
frequency  at  which  all  of  the  simulations  were  conducted  was  20  MHz  [38] .  The  inhomo¬ 
geneity  parameters  are  varied  according  to  Table  4.3.  Each  perturbation  case  includes  an 
absorption  coefficient  of  100,  200,  300,  or  400%  of  0.05  cm-1,  or  rather  a  contrast  of  0 
(none),  100  (low),  200  (high),  or  300%  (very  high)  relative  to  the  background  absorption 
coefficient,  respectively.  It  should  be  noted  that  these  absorption  contrasts  are  still  less 
than  actual  cancerous  tumor  values  given  in  Table  4.3.  Similarly,  the  scattering  parameter 
is  100,  110,  150,  or  200%  of  10  cm-1,  meaning  a  contrast  of  0  (none),  10  (low),  50  (high), 
or  100%  (very  high)  relative  to  the  background  reduced  scattering  coefficient.  A  pure 
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absorber  indicates  an  inhomogeneity  with  only  an  absorption  contrast,  while  a  pure  scat- 
terer  is  an  inhomogeneity  with  only  a  scattering  contrast.  From  the  simulated  results  of 
each  of  the  perturbation  cases  in  the  analytic  transfer  function  code  as  well  as  the  Fourier 
optics  method  code,  the  effects  the  different  levels  of  contrast  have  on  the  magnitude  and 
phase  change  in  the  resulting  total  wave  are  determined.  All  of  the  cases  in  Table  4.3  do 
not  exceed  the  conditions  for  the  Pi  approximation  of  the  linear  transport  equation,  e.g. 
Ma  <<  fj,'a  and  fia/[fis  +  Mo]  ~  1- 


Perturbation 

Absorption 

Scattering 

Mo, in 

Ms, in 

Case 

Contrast 

Contrast 

(cm-1) 

(cm-1) 

Base  Line 

none 

none 

0.05 

10 

Pure  Absorption 

high 

none 

0.10 

10 

Pure  Absorption 

high 

none 

0.15 

10 

Pure  Absorption 

very  high 

none 

0.20 

10 

Pure  Scatterer 

none 

high 

0.05 

11 

Pure  Scatterer 

none 

high 

0.05 

15 

Pure  Scatterer 

none 

very  high 

0.05 

20 

Mixed 

low 

low 

0.10 

11 

Mixed 

low 

high 

0.10 

15 

Mixed 

low 

very  high 

0.10 

20 

Mixed 

high 

low 

0.15 

11 

Mixed 

high 

high 

0.15 

15 

Mixed 

high 

very  high 

0.15 

20 

Mixed 

very  high 

low 

0.20 

11 

Mixed 

very  high 

high 

0.20 

15 

Mixed 

very  high 

very  high 

0.20 

20 

Table  4.3  The  wave  behavior  is  analyzed  for  highly  scattering  optical  parameters,  as 
well  as  for  levels  of  contrast  between  the  background  tissue  (designated  by  the 
subscript  “out”)  and  the  inhomogeneity  (designated  by  the  subscript  “in”). 
All  of  these  values  are  for  a  homogeneous  background  medium  with  fia,out  = 
0.05  cm-1  and  /4j0Ut  =  10  cm-1. 


4-10.2  Analysis  Plan.  The  behavior  of  a  weakly  perturbative  system  is  analyzed 
in  several  ways  using  the  cases  listed  in  Table  4.3.  In  addition,  the  extent  of  error  in  the 
Fourier  optics  approximation  is  also  studied  under  this  system. 

4-10.2.1  Sensitivity  Analysis.  For  the  low,  high,  and  mixed  contrast  cases, 
the  total  wave  behavior  relative  to  the  contrast  in  the  optical  characteristics  is  investigated. 
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The  moment  analysis  in  Chapter  III  determined  that  an  absorptive  contrast  results  in 
an  amplitude  change,  while  a  scattering  contrast  causes  perturbation  in  the  phase.  The 
sensitivity  analysis  will  show  to  what  extent  varying  degrees  of  contrast  perturb  the  output 
wave  amplitude  and  phase.  For  a  unit  amplitude,  zero  phase  incident  wave,  the  transfer 
function  in  Equation  (3.26)  can  be  used  to  examine  this  sensitivity  of  the  system. 

4-10.2.2  Size  of  Inhomogeneity.  A  weakly  perturbative  case  is  used  to  de¬ 
termine  the  effects  low  contrast  in  absorption  and  scattering  have  on  detecting  various  sizes 
of  the  inhomogeneity.  Again,  the  transfer  function  calculated  from  the  analytic  transfer 
function  model  can  be  used  in  this  analysis. 

4-10.2.3  Moment  Contribution.  For  each  perturbative  case,  the  moments  of 
the  scattered  wave  are  investigated  to  determine  the  amount  of  contribution  each  moment 
has  to  the  total  wave  relative  to  the  amount  of  absorptive  and/or  scattering  contrast 
present  in  the  system. 

4-10.2.4  Validity  of  Fourier  Optics  Approximation.  The  total  output  wave 
calculated  from  the  analytic  transfer  function  model  is  compared  to  Goodman’s  Fourier 
optics  method.  Since  the  Fourier  optics  approach  is  a  straight  ray  approximation  to  distor¬ 
tion  caused  by  a  spherical  inhomogeneity,  the  relative  error  in  the  two  models,  specifically 
in  the  model  of  the  inhomogeneity,  will  determine  the  regions  of  validity  of  this  approach. 

4-11  Summary 

Three  different  models  of  wave  propagation  through  a  homogeneous  medium  con¬ 
taining  a  spherical  object  were  simulated  in  this  chapter.  The  first  is  an  implementation 
of  the  analytic  solution  (derived  by  Boas)  to  the  problem.  The  second  is  a  Fourier  optics 
approach  that  applies  an  exact  formulation  of  a  transfer  function  to  propagate  the  wave 
from  the  source  to  the  detector.  The  last  is  a  plane  wave  representation  of  the  incident 
spherical  wave  developed  by  Goodman  and  uses  a  simpler  transfer  function  to  propagate 
the  perturbed  wave  from  the  inhomogeneity  to  the  detector. 
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The  first  two  models  were  correlated  against  a  known  working  code,  PMI.  With  the 
exception  of  a  complex  scale  factor  and  an  approximated  parameter,  the  models  are  in 
good  agreement  with  PMI.  Goodman’s  wave  propagation  model  is  examined  in  Chapter  V 
for  validity. 

The  plan  for  analysis  of  the  analytic  transfer  function  was  outlined  for  a  range  of 
weak  to  strong  perturbative  cases  in  Table  4.3.  The  surrounding  background  medium  is 
modeled  as  human  breast  tissue  with  a  fixed  set  of  optical  parameters.  By  varying  the 
optical  properties  of  the  inhomogeneity,  contrasts  between  the  properties  in  the  object 
and  the  surrounding  medium  create  the  weak  and  strong  perturbative  conditions.  In 
addition  to  varying  the  optical  properties,  changing  the  size  of  the  inhomogeneity  will  be 
examined.  By  measuring  the  amount  of  perturbation  caused  by  the  contrast  in  the  optical 
properties  or  by  the  size  of  the  object,  the  ability  to  detect  the  inhomogeneity  can  be 
determined.  The  amount  of  contrast  present  in  the  system  also  affects  which  moments 
significantly  contribute  to  the  perturbation  in  the  output  wave.  By  determining  which 
moments  dominate  in  each  perturbation  case,  an  improved  approximation  can  be  made  to 
the  transfer  function  expression  developed  in  Chapter  III. 
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V.  Simulation  Results  and  Analysis 


5.1  Introduction 

Tumors  are  detectable  because  their  optical  scattering  and  absorption  properties  are 
different  than  that  of  the  media  in  which  they  are  embedded.  This  contrast  in  optical 
properties  causes  a  perturbation  in  the  diffuse  photon  density  wave  (DPDW).  If  the  level 
of  change  in  the  DPDW  is  greater  than  a  certain  limit,  the  tumor  is  considered  detectable. 
An  understanding  of  the  interaction  and  the  levels  of  contrasts  needed  to  detect  and  localize 
the  object  is  important  optical  imaging,  particularly  in  medical  applications. 

General  behavior  of  the  total  wave  for  different  absorptive  and  scattering  proper¬ 
ties  of  the  inhomogeneity  relative  to  a  constant  background  medium  are  considered  in 
this  evaluation.  The  particular  environment  of  interest  in  this  thesis  is  a  weak  pertur¬ 
bation,  meaning,  the  optical  characteristics  of  the  inhomogeneity  slightly  differ  from  the 
background  medium.  Since  the  transfer  function  developed  in  Chapter  III  completely  char¬ 
acterizes  the  perturbations  caused  by  the  inhomogeneity,  the  perturbative  behavior  can  be 
determined  by  examining  the  transfer  function  in  the  spatial  domain.  The  contribution 
of  the  individual  moments  in  the  scattered  wave  are  also  studied  to  determine  the  neces¬ 
sary  number  of  moments  needed  to  best  approximate  the  scattered  wave.  This  chapter 
concludes  with  a  comparison  between  the  analytic  transfer  function  and  the  Fourier  optics 
model  in  Section  2.5  to  determine  the  regions  of  validity  as  well  as  the  relative  accuracy 
of  the  transmittance  function  in  Goodman’s  method. 

5.2  Sensitivity  to  Contrasts  in  Optical  Parameters 

Perturbations  caused  by  the  presence  of  an  inhomogeneity  are  caused  by  a  differ¬ 
ence,  or  contrast,  in  the  optical  parameters  inside  and  outside  the  object.  The  amount 
of  distortion  by  the  contrast  defines  the  sensitivity  of  the  system.  By  measuring  the  lev¬ 
els  of  distortion,  inhomogeneities  can  be  detected.  The  measurement  precision  necessary 
to  detect  optical  inhomogeneities  can  be  estimated  using  the  analytic  transform  solution 
method.  The  required  amplitude  precision  is  determined  by  the  ratio  of  the  magnitudes 
of  the  total  wave  to  the  homogeneous  wave.  The  difference  in  the  phase  between  the  total 
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Figure  5.1  The  transfer  function  in  the  spatial  domain  is  plotted  for  an  inhomogeneity 
with  n'$  =  15  cm-1  and  /ia  =  0.15  cm-1  embedded  in  a  medium  with  the 
n's  =  10  cm-1  and  (ia  =  0.05  cm-1.  This  shows  that  the  best  estimate 
for  measurement  precision  can  be  determined  here  at  the  x  =  0  position  for 
y  =  0  cm  and  the  detector  plane  located  at  z  =  3  cm. 

and  the  homogeneous  wave  is  the  necessary  phase  precision.  If  the  homogeneous  wave  has 
unit  amplitude  and  zero  phase,  the  magnitude  and  phase  precision  required  is  simply  the 
magnitude  and  phase  of  the  transfer  function  from  Equation  (3.26)  in  the  spatial  domain. 

5.2.1  System  Configuration.  To  study  the  system  sensitivity,  the  system  is  con¬ 
figured  as  in  Section  4.2.  Figure  5.1  illustrates  that  for  an  inhomogeneity  with  a  reduced 
scattering  coefficient,  fi's  =  15  cm-1,  and  an  absorption  coefficient,  fia  =  0.15  cm-1, 
embedded  in  a  medium  with  (i's  =  10  cm-1  and  ya  =  0.05  cm-1,  the  maximum  per¬ 
turbation  in  the  DPDW  occurs  at  the  x  =  0  position.  Since  the  system  is  symmetric,  the 
same  view  holds  for  the  y  direction.  The  locations  along  the  x  and  y  axes  indicate  that  the 
maximum  perturbation  occurs  along  the  axis  that  contains  the  source,  the  inhomogeneity 
and  the  detector.  Along  that  axis,  then,  the  best  estimate  for  measurement  precision  can 
be  determined. 

For  the  perturbation  cases  listed  in  Table  4.3,  the  precision  in  the  magnitudes  and 
phases  are  determined  and  compared  to  what  most  detectors  can  detect,  which  is  at  least 
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a  0.1%  amplitude  and  0.1°  phase  change  [10].  If  the  perturbation  caused  by  the  inhomo¬ 
geneity  is  greater  than  the  detector  level,  the  object  is  considered  detectable.  It  should 
be  noted  that  noise  is  not  added  to  the  signal  and  that  detection  ultimately  depends  on 
the  level  of  noise  in  the  system  [13].  For  this  reason,  the  amplitude  detection  level  con¬ 
sidered  in  this  analysis  will  be  increased  to  1.0%.  For  all  the  contrast  cases,  the  system 
is  configured  as  in  Section  4.2  at  a  frequency  of  20  MHz.  The  spherical  inhomogeneity  is 
embedded  in  a  background  medium  with  p's  =  10  cm-1  and  pa  —  0.05  cm-1,  and  the 
optical  parameters  of  the  inhomogeneity  are  varied. 

In  addition  to  being  able  to  detect  the  inhomogeneity,  the  minimum  in  either  the 
magnitude  or  the  phase  of  the  transfer  function  in  the  spatial  domain  can  be  used  to 
determine  the  depth  of  the  object. 

5.2.2  Absorption  Contrast  With  a  Pure  Absorber.  A  pure  absorber  with  p's  = 
10  cm-1  is  studied  in  this  analysis.  The  absorption  coefficient  for  the  absorber  is  varied  to 
show  the  perturbation  effects  on  the  DPDW. 

In  Figure  5.2,  a  slice  along  the  (x,  y)  —  (0, 0)  is  shown  for  the  magnitude  and  phase  of 
the  homogeneous  wave  through  the  media,  the  scattered  wave  caused  by  the  absorber,  and 
the  total  wave  in  the  system.  It  is  along  this  central  axis  that  the  maximum  perturbation 
will  occur  as  discussed  in  Section  5.2.1. 

The  plot  in  Figure  5.3  indicates  for  a  pure  absorber  at  a  detector  plane  3  cm  from 
the  center  of  the  inhomogeneity,  as  small  as  4%  amplitude  precision  is  necessary  to  detect 
an  inhomogeneity  with  100%  absorption  contrast.  There  is  no  detectable  phase  difference 
at  that  distance.  Note  that  the  straight  solid  line  indicates  the  detection  level,  but  the 
discontinuous  solid  line  is  for  an  inhomogeneity  with  no  contrast  in  either  scattering  or 
absorption,  and  consequently  no  noticeable  change  in  phase  or  magnitude.  As  the  value 
of  the  absorption  coefficient  increases  for  the  inhomogeneity  relative  to  the  background 
medium,  the  precision  required  is  less  stringent.  For  example,  at  three  times  the  amount 
of  absorption  contrast,  only  8%  amplitude  precision  is  needed. 

If  the  detector  plane  moves  closer  to  the  inhomogeneity,  the  relative  amount  increase 
in  both  magnitude  and  phase  are  evident.  However,  even  at  the  boundary  of  the  inhomo- 


5-3 


Figure  5.2  A  view  along  the  ( x ,  y)  =  (0, 0)  axis  through  the  medium  is  shown  for  a  pure 
absorber  with  an  absorption  coefficient,  fia  =  0.15  cm-1.  The  (a)  magnitude 
and  (b)  phase  of  the  waves  shown  are  the  homogeneous  (dashed  line),  the 
scattered  wave  due  to  the  inhomogeneity  (dashed-dotted  line),  and  the  total 
wave  (solid  line).  It  is  along  this  axis  that  the  maximum  perturbation  will 
occur. 
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geneity,  the  phase  is  only  0.08°  and  does  not  meet  the  minimum  detection  level  for  most 
detectors. 

5.2.3  Scattering  Contrast  With  a  Pure  Scatterer.  In  this  case,  a  pure  scatterer 
with  na  =  0.05  cm-1  is  investigated.  To  show  the  perturbation  effects  on  the  DPDW, 
the  reduced  scattering  coefficient  for  the  object  is  varied  from  no  scattering  contrast  up  to 
twice  the  value  of  the  background  medium. 

In  Figure  5.4,  a  slice  along  the  ( x ,  y)  —  (0, 0)  is  shown  for  the  magnitude  and  phase 
of  the  homogeneous  wave  through  the  medium,  the  scattered  wave  caused  by  the  scatterer, 
and  the  total  wave  in  the  system.  As  discussed  in  Section  5.2.1,  the  maximum  perturbation 
will  occur  along  this  axis. 

The  plot  in  Figure  5.5  shows  that  at  the  detector  plane,  3%  amplitude  and  0.03° 
phase  precision  are  necessary  to  detect  the  inhomogeneity  at  a  3  cm  distance  from  the 
object  center,  with  a  10%  scattering  contrast  present.  The  phase  amount  is  not  sufficient 
to  meet  most  detector  minimum  levels.  If  the  detector  plane  is  located  close  to  the  right 
boundary  of  the  inhomogeneity,  the  phase  perturbation  increases  to  0.1°.  The  maximum 
magnitude  and  phase  in  all  variations  of  the  absorption  contrast  occur  at  the  right  (positive 
z )  boundary  of  the  inhomogeneity. 

As  the  scattering  contrast  increases  for  the  inhomogeneity  relative  to  the  background 
medium,  the  perturbation  in  both  the  magnitude  and  phase  increases.  At  a  scattering 
contrast  of  1.5  times  the  background  scattering  value,  the  phase  is  0.15°  at  the  detector 
plane  3  cm  from  the  center  of  the  inhomogeneity  and  meets  the  lower  detector  limits.  Note 
that  the  straight  solid  line  indicates  the  detection  level,  but  the  discontinuous  solid  line  in 
Figure  5.5  is  for  an  inhomogeneity  with  no  contrast  in  either  scattering  or  absorption,  and 
consequently  no  noticeable  phase  or  magnitude. 

The  results  from  Sections  5.2.2  and  5.2.3  are  consistant  with  results  presented  in  the 
literature  for  purely  absorbtive  and  purely  scattering  inhomogeneities  [12,13]. 

5.2.4  Low  Contrast  in  Absorption  and  Varied  Contrast  in  Scattering  Coefficients. 

In  this  case,  an  inhomogeneity  with  a  100%  absorption  contrast,  or  pa,in  —  0.10  cm-1,  is 
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Figure  5.3  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  a  pure  absorber  along  the  z-&xis. 
The  object  absorption  coefficient  is  varied  as:  /i0)jn  =  0.05  cm-1  (solid  line), 
Ha, in  =  0-10  cm-1  (dashed  line),  \ia^n  =  0.15  cm-1  (dashed-dotted  line), 
and  [Aatin  =  0.20  cm-1  (dotted  line).  For  both  figures,  fi'sin  =  10  cm-1, 
Ha, out  —  0.05  cm-1,  and  n's,out  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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Figure  5.4  A  view  along  the  (x,  y )  =  (0, 0)  axis  through  the  medium  is  shown  for  a  pure 
scatterer  with  reduced  scattering  coefficient,  y!s  =  15  cm-1.  The  waves 
shown  are  the  homogeneous  (dashed  line),  the  scattered  wave  due  to  the 
inhomogeneity  (dashed-dotted  line),  and  the  total  wave  (solid  line).  It  is 
along  this  axis  that  the  maximum  perturbation  will  occur. 
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Figure  5.5  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  a  pure  scatterer  along  the  z-axis. 
The  object  scattering  coefficient  is  varied  as:  n's  in  =  10  cm-1  (solid  line), 
Ms, in  =  11  cm-1  (dashed  line),  n's  in  =  15  cm-1  (dashed-dotted  line), 
and  n's>in  =  20  cm-1  (dotted  line).  For  both  figures,  najn  =  0.05  cm-1, 
Ma,out  =  0-05  cm-1,  and  (i'Si0Ut  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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studied.  The  reduced  scattering  coefficient  for  the  object  is  varied  to  examine  the  resulting 
perturbations  for  this  low  absorption  contrast.  Figure  5.6  indicates  that  for  a  10%  difference 
in  scattering  coefficients,  4%  amplitude  and  0.04°  phase  precision  are  necessary  to  detect 
the  inhomogeneity  at  the  detector  plane  3  cm  from  the  center  of  the  inhomogeneity.  At 
this  low  level  of  scattering  contrast,  there  is  insufficient  phase  perturbation  to  be  detected. 

As  in  the  pure  scattering  case,  an  increase  in  the  scattering  contrast  decreases  the 
amount  of  precision  required  to  detect  the  object.  Note  that  the  straight  solid  line  indicates 
the  detection  level,  but  the  discontinuous  solid  line  in  the  graph  illustrates  no  scattering 
contrast.  The  absorption  contrast  causes  a  magnitude  perturbation  as  seen  before  in 
the  pure  absorber  case.  When  a  scattering  contrast  is  present,  the  resulting  magnitude 
perturbation  is  in  addition  to  the  perturbation  caused  by  the  absorption  contrast.  A 
similar  effect  is  seen  in  the  phase.  However,  since  the  phase  due  to  the  absorption  contrast 
is  negative,  but  is  positive  due  to  the  scattering  contrast,  the  combination  of  the  two  is 
a  lower  phase  amount  than  if  there  was  no  absorption  contrast  present.  In  other  words, 
the  phase  perturbation  due  to  an  absorption  contrast  combines  destructively  with  the 
perturbation  due  to  a  scattering  contrast. 

The  maximum  in  both  the  magnitude  and  phase  occurs  at  the  inhomogeneity  bound¬ 
ary  at  0.5  cm.  There  is  an  additional  0.2°  phase  drop-off  to  the  right  of  the  object  boundary 
due  to  the  contrast  in  absorption. 

5.2.5  High  Contrast  in  Absorption  and  Varied  Contrast  in  Scattering  Coefficients. 

The  inhomogeneity  absorption  coefficient  is  increased  to  three  times  the  background 
medium  value,  meaning  //0)in  =  0.15  cm-1.  Again,  the  reduced  scattering  coefficient  for  the 
scatterer  is  varied  to  show  the  perturbations  for  high  absorption  contrast  in  the  resulting 
wave.  The  plot  in  Figure  5.7  illustrates  that  at  the  detector  plane  at  z  =  3  cm,  as  small 
as  8%  amplitude  and  0.04°  phase  precision  is  necessary  to  detect  the  inhomogeneity  with 
10%  scattering  contrast.  Moving  the  detector  plane  does  not  increase  the  phase  to  a 
detectable  level  at  this  scattering  contrast  amount.  If  the  scattering  contrast  is  increased 
to  1.5  times  that  of  the  background,  the  phase  precision  increases  to  0.12°  and  is  sufficient 
to  be  detected.  The  corresponding  magnitude  precision  increases  only  to  11%.  As  in 
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Figure  5.6  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  ana¬ 
lytic  transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  low 
contrast  in  absorption  and  varied  contrast  in  scattering  along  the  2-axis. 
The  object  scattering  coefficient  is  varied  as:  ns^n  =  10  cm-1  (solid  line), 
^'s,in  =  cm-1  (dashed  line),  n'ain  =  15  cm-1  (dashed-dotted  line), 
and  fi'a  in  =  20  cm-1  (dotted  line).  For  both  figures,  jia,in  =  0.10  cm-1, 
Ha, out  —  0-05  cm-1,  and  fi's,out  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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the  previous  case  of  a  small  amount  of  absorption  contrast,  the  scattering  and  absorption 
contrast  effects  combine  to  cause  a  higher  detectable  amplitude  and  a  lower  detectable 
phase. 

The  location  of  the  maximum  in  the  phase  occurs  at  the  inhomogeneity  boundary  at 
0.5  cm.  Any  additional  phase  drop-off  to  the  right  of  the  boundary  is  not  noticeable. 

Figure  5.8  shows  the  sensitivity  of  an  inhomogeneity  modeled  with  a  very  high  ab¬ 
sorption  coefficient  of  payin  =  0.20  cm-1,  or  four  times  that  of  the  background  medium. 
For  a  10%  contrast  in  scattering,  the  necessary  magnitude  precision  is  12%  ,  and  phase 
precision  is  0.04°.  Only  the  magnitude  perturbation  will  be  detected  by  most  detectors. 

Moving  the  detector  plane  towards  the  inhomogeneity,  both  the  magnitude  and  phase 
increase.  Yet,  even  at  the  boundary  of  the  inhomogeneity  at  10%  contrast  in  scattering,  the 
phase  remains  undetectable.  However,  scattering  contrasts  of  50%  and  100%  are  detectable 
in  the  entire  region  between  the  inhomogeneity  boundary  and  the  detector  plane  at  z=3 
cm.  Consistent  with  the  previous  absorption  contrast  cases,  the  scattering  and  absorption 
contrasts  effects  combine  to  cause  a  higher  detectable  amount  in  the  amplitude  and  a 
lower  detectable  amount  in  the  phase.  However  the  negative  phase  contribution  due  to  the 
absorption  contrast  is  minimal  when  compared  to  the  phase  amount  due  to  the  scattering 
contrast. 

5.2.6  Low  Contrast  in  Scattering  and  Varied  Contrast  in  Absorption  Coefficients. 
The  sensitivity  of  an  inhomogeneity  with  a  10%  increase  in  scattering,  or  p!s  in  =  11  cm-1, 
is  investigated.  To  show  the  resulting  perturbation  in  the  DPDW,  the  absorption  coefficient 
for  the  object  is  varied  from  no  contrast  to  four  times  the  value  of  the  background  medium. 
The  plot  in  Figure  5.9  indicates  that  for  an  absorption  contrast  of  100%,  5%  amplitude 
and  0.03°  phase  precision  is  necessary  to  detect  the  inhomogeneity  at  the  detector  plane 
3  cm  from  the  center  of  the  inhomogeneity.  The  phase  is  not  sufficient  to  be  detected 
by  most  detectors.  However,  the  amplitude  is  well  above  the  necessary  limits.  As  in 
the  pure  absorption  case,  an  increase  in  the  absorption  contrast  decreases  the  amount  of 
magnitude  precision  required  to  detect  the  object.  With  the  small  amount  of  scattering 
contrast  present,  the  resulting  magnitude  perturbation  is  in  addition  to  that  caused  by 
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Figure  5.7  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  ana¬ 
lytic  transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  high 
contrast  in  absorption  and  varied  contrast  in  scattering  along  the  i-axis. 
The  object  scattering  coefficient  is  varied  as:  =  10  cm-1  (solid  line), 

fj,'s  in  =  11  cm-1  (dashed  line),  n's  in  =  15  cm-1  (dashed-dotted  line), 

and/x'sin  =  20  cm-1  (dotted  line).  For  both  figures,  iia>in  =  0.15  cm-1, 
fia,out  =  0-05  cm-1,  and  ^'Sj0ut  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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Figure  5.8 


These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  very  high 
contrast  in  absorption  and  varied  contrast  in  scattering  along  the  2-axis. 
The  object  scattering  coefficient  is  varied  as:  \)!s  in  =  10  cm-1  (solid  line), 
/z')in  =  11  cm-1  (dashed  line),  n'ain  =  15  cm-1  (dashed-dotted  line), 
and  =  20  cm-1  (dotted  line).  For  both  figures,  fia,in  =  0.20  cm-1, 

Ha, out  =  0.05  cm-1,  and  pJs>out  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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increased  absorption  contrast.  A  similar  effect  is  seen  in  the  phase  on  the  right  side 
of  the  inhomogeneity  (i.e.  along  the  positive  z  direction).  The  magnitude  perturbation  is 
detectable  for  all  levels  of  absorption  contrast  regardless  of  the  detector  position.  However, 
the  corresponding  phase  perturbation  is  not  detectable  at  any  point  in  the  region. 

5.2.7  High  Contrast  in  Scattering  and  Varied  Contrast  in  Absorption  Coefficients. 

The  scattering  contrast  is  increased  to  50%,  meaning  the  inhomogeneity  has  a  scattering 

coefficient  of  p's  in  =  15  cm-1.  Again,  the  absorption  coefficient  for  the  object  is  varied 
to  show  the  perturbation  effects  on  the  DPDW.  Figure  5.10  illustrates  that  in  the  region 
between  the  inhomogeneity  boundary  and  the  boundary  of  the  system  at  3  cm,  both  the 
magnitude  and  the  phase  amounts  are  detectable  for  all  absorption  contrast  amounts. 
The  precision  required  to  detect  a  change  in  the  magnitude  or  the  phase  decreases  as  the 
detector  plane  is  moved  towards  the  center  of  the  inhomogeneity.  Additional  increase  in 
absorption  contrast  does  not  significantly  perturb  the  phase  further. 

Figure  5.11  shows  the  sensitivity  of  an  inhomogeneity  modeled  with  a  very  high 
scattering  coefficient  of  p,'s  in  =  20  cm-1,  or  double  that  of  the  background  medium.  As  in 
the  previous  scattering  contrast  case,  both  the  magnitude  and  phase  amounts  are  detectable 
for  all  absorption  contrast  variations  in  the  region  bounded  by  the  inhomogeneity  and  the 
detector  plane  at  z=3  cm. 

Consistent  with  the  previous  scattering  contrast  cases,  the  scattering  and  absorption 
contrasts  effects  combine  to  cause  a  higher  detectable  amount  in  the  amplitude.  However, 
the  absorption  coefficient  has  little  effect  on  the  detectability  of  the  phase. 

5.2.8  Summary  of  Sensitivity  Analysis.  In  this  section,  the  amount  of  perturba¬ 
tion  in  the  total  wave  was  examined  for  various  combinations  of  optical  contrasts  between 
the  inhomogeneity  and  the  surrounding  medium.  The  results  show  that  a  high  absorbtive 
contrast  causes  a  decrease  in  the  magnitude  of  the  perturbed  wave  but  does  not  signifi¬ 
cantly  affect  the  phase.  However,  a  scattering  contrast  causes  a  significant  change  in  the 
phase  of  the  perturbed  wave  but  minimally  affects  the  magnitude.  A  combination  of  two 
contrasts  combine  their  effects  in  the  magnitude  and  the  phase  of  the  total  wave. 
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Figure  5.9  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  low  con¬ 
trast  in  scattering  and  varied  contrast  in  absorption  along  the  z-axis.  The 
object  absorption  coefficient  is  varied  as:  na^n  =  0.05  cm-1  (solid  line), 
Mo, in  =  0.10  cm-1  (dashed  line),  \ia^n  =  0.15  cm-1  (dashed-dotted  line), 
and  Mo, in  =  0.20  cm-1  (dotted  line).  For  both  figures,  ii's  in  =  11  cm-1, 
Mo, out  =  0.05  cm-1,  and  =  10  cm-1.  The  straight  solid  line  in  each 

graph  indicates  the  minimum  detection  level. 
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Figure  5.10  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  high  con¬ 
trast  in  scattering  and  varied  contrast  in  absorption  coefficients  along  the 
z-axis.  The  object  absorption  coefficient  is  varied  as:  =  0.05  cm-1 

(solid  line),  /x0)jn  =  0.10  cm-1  (dashed  line),  na^n  —  0.15  cm-1  (dashed- 
dotted  line),  and  na^n  —  0.20  cm-1  (dotted  line).  For  both  figures, 

Vs,in  =  15  cm_1)  Va,out  =  0.05  cm-1,  and  /z')0ttt  =  10  cm-1.  The 

straight  solid  line  in  each  graph  indicates  the  minimum  detection  level. 
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(b) 

Figure  5.11  These  plots  show  the  (a)  magnitude  and  (b)  phase  (degrees)  of  the  analytic 
transfer  function  in  the  spatial  domain  for  an  inhomogeneity  with  very  high 
contrast  in  scattering  and  varied  contrast  in  absorption  along  the  2-axis.  The 
object  absorption  coefficient  is  varied  as:  iia>in  =  0.05  cm-1  (solid  line), 
Mo.in  =  0.10  cm-1  (dashed  line),  [ia.in  =  0.15  cm-1  (dashed-dotted  line), 
and  na)in  =  0.20  cm-1  (dotted  line).  For  both  figures,  ii'sin  =  20  cm-1, 
Va,ov,t  =  0.05  cm-1,  and  n'S)OUt  =  10  cm-1.  The  straight  solid  line  in  each 
graph  indicates  the  minimum  detection  level. 
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The  sensitivity  results  also  showed  that  the  location  of  the  detector  plane  relative 
to  the  inhomogeneity  determines  the  amount  of  perturbation  detected.  The  maximum 
change  in  both  the  amplitude  and  the  phase  in  the  total  wave  occured  at  the  boundary  of 
the  inhomogeneity,  assuming  the  object  is  between  the  source  and  detector.  This  maximum 
position  can  be  used  to  determine  the  depth  of  the  object. 

Differing  optical  parameters  are  only  one  perturbative  effect.  In  the  next  section  the 
disruption  in  the  DPDW  caused  by  variations  in  the  size  of  the  inhomogeneity  is  examined. 
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5.3  Sensitivity  to  Inhomogeneity  Size 

The  perturbation  of  the  incident  wave  due  to  an  inhomogeneity  is  directly  depen¬ 
dent  on  the  size  of  the  inhomogeneity.  The  analytic  analysis  in  Section  2.4.1  showed  this 
relationship.  It  has  been  shown  that  for  an  absorption  contrast  of  three  times  that  of 
the  background  medium,  or  a  scattering  contrast  factor  of  1.5,  an  inhomogeneity  can  be 
detected  [13].  For  a  mix  of  low  contrast  in  absorption  and  scattering  relative  to  the  back¬ 
ground  medium,  size  is  a  critical  factor  in  perturbing  the  wave  to  a  level  of  detection 
ability.  The  system  is  configured  as  described  in  Section  5.2.1.  Figure  5.12  illustrates  the 
sensitivity  levels  of  different  sizes  of  inhomogeneities  that  have  weakly  perturbing  optical 
characteristics  of  100%  in  absorption  contrast  and  10%  in  scattering  contrast  with  the 
background  medium,  namely  y,a  =  0.10  cm-1  and  y!s  =  11  cm-1. 

At  the  detector  plane,  for  a  detector  minimum  limit  of  1.0%  magnitude  perturbation 
(indicated  by  the  straight  solid  line),  only  the  6  mm,  8  mm  and  10  mm  diameter  objects 
can  be  decisively  detected.  For  a  minimum  limit  of  0.1°  phase  (also  indicated  by  the 
straight  solid  line),  none  of  the  objects  are  detectable.  If  the  detector  plane  is  moved  to 
the  peak  amount  of  each  object,  all  save  the  1  mm  diameter  object  can  be  detected  by 
magnitude.  Only  the  10  mm  object  is  detected  on  phase  perturbation  alone  and  requires 
the  detector  to  be  on  the  same  side  of  the  object  as  the  source,  consequently  violating  the 
system  geometry. 

5-4  Scattering  Moment  Contributions 

In  Chapter  II,  the  scattering  wave  was  shown  to  be  comprised  of  a  series  of  moments 
that  depended  on,  among  other  quantities,  the  contrast  of  optical  properties.  If  the  condi¬ 
tions  \kouta\  <<  1  and  \kina\  «  1  are  valid,  where  kout  (fcjn)  is  the  wavenumber  outside 
(inside)  the  inhomogeneity  with  radius  a,  only  the  first  three  moments  of  the  series  are 
detectable  in  experiments  due  to  signal-to- noise  limitations  [13]. 

However,  if  the  conditions  mentioned  above  are  not  met,  additional  moments  are 
detectable.  For  the  contrast  cases  listed  in  Table  4.3  and  for  a  1  cm  diameter  inhomogeneity, 
\kina\  =  1.069,  and  \kouta\  ranges  from  0.6143  to  0.6158.  These  values  do  not  conclusively 
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1.05 


(b) 

Figure  5.12  These  sensitivity  plots  show  the  (a)  magnitudes  and  (b)  phases  (degrees) 
of  the  transfer  function  in  the  spatial  domain  with  low  contrast  in  scat¬ 
tering  absorption  coefficients  along  the  5-axis.  The  diameter  of  the  inho¬ 
mogeneity  is  varied  as:  10  mm  (solid  line),  8  mm  (dashed  line),  6  mm 
(dashed-dotted  line),  and  4  mm  (dotted  line),  2  mm  (>),  1  mm  (+).  For 
both  figures,  /xa, in  =  0.10  cm-1,  fi's  in  =  11  cm-1,  na>out  =  0.05  cm-1, 
and  pi'Sj0Ut  =  10  cm-1.  The  straight  solid  line  in  each  graph  indicates  the 
minimum  detection  level. 
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meet  the  conditions,  so  scattered  wave  moments  in  addition  to  the  first  three  are  expected 
to  have  significant  contributions  to  the  total  perturbation  caused  by  the  inhomogeneity. 

The  first  moment  (monopole)  was  shown  in  Section  2.4.1  to  be  predominantly  affected 
by  a  change  in  absorptive  contrasts,  while  the  second  (dipole)  and  third  (quadrupole) 
moments  were  affected  by  scattering  contrasts.  In  Section  5.2.1,  the  sensitivity  analysis 
determined  that  although  contrasts  in  either  absorption  or  scattering  caused  a  detectable 
amount  in  the  amplitude,  the  absorption  contrast  caused  a  higher  detectable  amount  than 
the  scattering  contrast.  In  addition,  the  sensitivity  results  are  consistent  in  the  moment 
analysis  of  a  pure  absorber  and  pure  scatterer,  as  well  as  mixed  contrast  cases. 

In  this  section,  the  contributions  of  the  first  six  moments  in  the  scattered  wave  are 
investigated  for  several  of  the  perturbative  cases  listed  in  Table  4.3.  The  contribution  of  an 
individual  moment  is  determined  by  the  absolute  value  of  the  ratio  of  the  magnitude  and 
phase  of  the  scattered  wave  moment  and  the  incident  wave.  The  particular  cases  considered 
in  this  section  include  a  pure  absorber,  a  pure  scatterer,  and  inhomogeneities  with  low  and 
high  absorptive  and  low  and  high  scattering  contrasts.  The  very  high  contrast  cases  are 
not  investigated  since  the  focus  of  this  research  is  on  weak  perturbations  caused  by  low 
contrasts  in  the  optical  properties. 

5.4-1  Pure  Absorption.  In  Figures  5.13  and  5.14,  a  pure  absorber  with  a  high 
and  low  absorption  contrast  are  illustrated,  respectively.  Overall,  a  higher  contrast  in 
absorption  causes  a  greater  amount  of  perturbation  in  the  magnitude  of  the  wave  than  a 
lower  contrast  case  causes.  Specifically,  for  an  inhomogeneity  with  an  absorption  coefficient 
three  times  that  of  the  background  medium  induces  twice  the  amount  of  amplitude  change 
than  an  inhomogeneity  with  an  absorption  coefficient  only  twice  that  of  the  background 
medium.  The  dominant  moment  for  a  pure  absorber  is  the  monopole.  Each  additional 
moment  contributes  to  the  magnitude  perturbation  on  the  order  of  a  factor  of  10  less  per 
moment.  For  example,  the  dipole  moment  is  a  factor  of  10  less  than  the  monopole  but  the 
quadrupole  is  a  factor  of  100  less  than  the  monopole. 
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As  expected  from  the  sensitivity  analysis  of  Section  5.2.2,  various  levels  of  absorp¬ 
tion  contrast  have  negligible  affect  on  the  phase  perturbations.  The  contributions  of  the 
moments  in  order  of  the  most  significant  is  the  first  through  the  sixth  moments. 

5-4-2  Pure  Scatterer.  For  a  pure  scatterer,  there  are  similar  trends  in  amplitude 
change  due  to  increased  contrast,  as  in  the  pure  absorption  case.  Figures  5.15  and  5.16 
depict  high  and  low  contrast  cases  for  a  pure  scatterer,  respectively.  The  high  contrast 
case  describes  an  inhomogeneity  with  a  reduced  scattering  coefficient  50%  higher  than  the 
background  medium.  The  high  contrast  causes  an  average  increase  of  3.8  times  the  amount 
of  amplitude  change  that  the  low  contrast  case  causes,  where  a  low  scattering  contrast  is 
defined  by  an  inhomogeneity  with  a  reduced  scattering  coefficient  only  10%  higher  than 
the  surrounding  medium.  The  contributions  by  moments  in  the  pure  scatterer  case  are 
significantly  different  than  in  the  pure  absorptive  case.  If  the  detector  plane  is  more  than 
a  radius  distance  from  the  inhomogeneity  boundary,  the  second  and  third  moments  are 
the  major  contributors  to  the  perturbation.  The  first  moment  ranks  third.  However, 
consistent  in  both  the  low  and  high  scattering  contrast  cases,  within  a  radius  distance  of 
the  inhomogeneity,  the  fourth  and  fifth  moments  contribute  more  than  the  monopole.  This 
means  that  in  the  detection  of  pure  scatterers,  the  position  of  the  detection  plane  is  an 
important  consideration  when  determining  how  many  moments  to  include  in  the  infinite 
series  representing  the  scattered  wave  (Equation  (2.11)). 

There  is  no  noticeable  phase  change  in  the  contributions  from  the  individual  moments 
due  to  an  increase  in  scattering  contrast.  However,  the  monopole  contributes  the  least 
amount  to  the  phase  perturbation,  as  expected  from  the  sensitivity  results  in  Section  5.2.3. 
An  interesting  behavior  of  the  monopole  shows  that  the  amount  of  its  phase  contribution 
depends  on  the  detector  plane  location. 

5-4-3  High  Absorption  Contrast.  The  presence  of  a  high  absorption  and  two 
different  levels  (high  and  low)  of  scattering  contrast  cause  different  amounts  of  change 
in  the  magnitude  of  the  perturbed  wave.  The  contributions  of  the  moments  also  vary 
depending  on  the  amount  of  absorptive  contrast  relative  to  the  scattering  contrast  present 
in  the  system.  Figures  5.17  and  5.18  illustrate  these  mixed  contrast  cases. 
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Figure  5.17  depicts  a  high  absorptive,  high  scattering  contrast  case.  Respectively,  the 
inhomogeneity  coefficients  are  a  factor  of  3  and  a  factor  of  1.5  of  the  coefficient  values  of 
the  surrounding  medium.  The  first  moment  is  the  major  contributer  to  the  perturbation, 
while  the  contributions  of  the  next  five  moments  monotonically  decrease.  The  presence 
of  the  high  scattering  contrast  relative  to  no  contrast  (Figure  5.13)  causes  an  increase  in 
contributions  from  the  third  through  the  sixth  moments  by  average  factors  of  9,  30,  50,  and 
150,  respectively.  Note  that  near  the  inhomogeneity,  the  second  moment  is  not  affected 
by  the  increase  in  the  scattering  contrast.  This  contribution  variation  indicates  the  second 
moment  depends  on  the  detector  plane  location. 

Keeping  a  high  absorptive  contrast  present  but  lowering  the  scattering  contrast  so 
the  inhomogeneity  scattering  coefficient  is  only  a  factor  of  1.1  of  the  coefficient  of  the 
surrounding  medium,  the  resulting  moment  contributions  to  the  magnitude  perturbation 
are  depicted  in  Figure  5.18.  The  respective  contributions  of  the  second  through  sixth 
moments  are  not  as  great  as  that  in  the  high  scattering  contrast  case.  However,  the 
contribution  of  the  second  moment  is  dependent  on  the  detector  plane  location.  Overall, 
the  first  moment  dominates  the  magnitude  perturbation.  At  a  minimum,  the  monopole  is 
26  times  the  amount  of  the  dipole,  which  is  the  next  highest  contributing  moment. 

Figures  5.17  and  5.18  depicting  an  inhomogeneity  with  a  high  absorptive  contrast 
relative  to  the  background  medium,  show  that  increasing  the  scattering  contrast  causes 
a  corresponding  increase  in  the  contributions  from  the  third  through  the  sixth  moments. 
The  second  moment  contribution  increases  when  the  detector  plane  is  not  close  to  the 
boundary  of  the  inhomogeneity.  Overall,  the  monopole  is  the  major  contributor  to  the 
perturbation  in  the  amplitude. 

The  presence  of  absorption  causes  the  monopole  to  also  be  the  major  contributor  to 
the  phase  perturbations.  However,  increasing  the  inhomogeneity  scattering  coefficient  from 
110%  to  150%  of  the  background  medium  (low  to  high  contrast)  causes  a  0.05°  decrease 
in  the  second  moment  and  a  0.01°  increase  in  the  third  moment.  The  remaining  moments 
are  unaffected  by  the  increase  in  the  scattering  contrast. 
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5.4-4  Low  Absorptive  Contrast.  Similar  to  the  cases  involving  a  high  absorptive 
contrast,  magnitude  perturbations  are  additionally  affected  by  increasing  scattering  con¬ 
trast  levels.  The  low  absorptive  contrast  with  low  and  high  scattering  contrast  cases  are 
graphed  in  Figures  5.19  and  5.20,  respectively. 

The  dominant  contributor  to  the  magnitude  perturbations  is  the  first  moment.  How¬ 
ever,  a  low  scattering  contrast  as  shown  in  Figure  5.19  causes  an  average  decrease  of  a  factor 
of  10  in  the  contribution  of  the  second  moment  as  compared  to  the  case  with  no  scattering 
present  (Figure  5.14).  Increasing  the  inhomogeneity  scattering  coefficient  to  150%  that  of 
the  surrounding  medium,  the  contributions  from  the  second  through  sixth  moments  are 
on  the  order  of  20  times  higher  than  in  the  no  scattering  case  as  seen  in  Figure  5.20.  Low 
scattering  contrast  in  the  presence  of  a  low  absorption  contrast  significantly  decreases  the 
dipole  moment,  but  that  effect  is  overcome  by  increasing  the  scattering  contrast. 

The  monopole  also  is  the  major  contributor  to  the  phase  perturbations  due  to  the 
presence  of  an  absorption  contrast.  An  increase  in  the  scattering  contrast  causes  an  increase 
in  the  second  moment  by  0.05°.  The  remaining  moments  are  unaffected. 

5.4.5  High  Scattering  Contrast.  Consider  the  cases  with  a  fixed  high  scatter¬ 
ing  contrast.  The  effects  caused  by  increasing  the  absorption  contrast  are  graphed  in 
Figures  5.20  and  5.17,  respectively. 

In  a  pure  scattering  case  as  shown  before  in  Figure  5.15,  the  dipole  moment  is  the 
major  contributor  to  the  magnitude  perturbations.  When  a  low  absorption  contrast  is 
introduced  into  the  system,  the  monopole  becomes  the  dominant  moment  (see  Figure  5.20). 
The  second  through  sixth  moments  are  not  affected  by  a  low  absorption  contrast  of  100% 
that  of  the  surrounding  medium.  As  the  absorption  contrast  increases  to  a  level  of  200% 
that  of  the  background  medium,  a  corresponding  increase  by  a  factor  of  2  in  the  monopole 
moment  is  apparent  in  Figure  5.17.  The  dipole  moment  decreases  by  an  average  factor  of 
0.8,  and  moments  three  through  six  are  not  noticeably  changed. 

In  the  presence  of  high  scattering  contrast,  an  increase  in  absorption  contrast  signif¬ 
icantly  increases  the  magnitude  perturbations  contributed  by  the  first  moment,  and  to  a 
lesser  degree,  decreases  the  third  moment. 
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Increasing  the  inhomogeneity  absorption  coefficient  from  100%  to  200%  that  of  the 
background  medium  (low  to  high  contrast)  does  not  affect  the  phase  perturbation  moment 
contributions.  This  behavior  is  expected  since  the  results  from  Section  5.2.2  showed  that 
in  the  presence  of  a  high  scattering  contrast,  the  absorption  contrast  does  not  significantly 
affect  the  distortion  in  the  phase  of  the  perturbed  wave. 

5.4.6  Low  Scattering  Contrast.  Decreasing  the  reduced  scattering  coefficient  of 
the  inhomogeneity  to  a  10%  contrast  with  the  background  medium  increases  the  effect 
that  an  absorption  contrast  has  on  the  magnitude  perturbation  caused  by  the  presence  of 
an  inhomogeneity.  Figures  5.18  and  5.19  illustrate,  respectively,  the  effects  high  and  low 
absorptive  contrasts  have  with  a  low  scattering  contrast. 

Similar  to  the  high  scattering  contrast  case  in  the  previous  section,  the  presence  of 
an  absorptive  contrast  causes  the  monopole  to  become  the  dominant  moment.  Comparing 
Figure  5.18  to  Figure  5.19  shows  that  increasing  the  absorption  contrast  200%  causes  a 
corresponding  doubling  of  the  amount  the  monopole  moment  contributes  to  the  magnitude 
perturbations.  The  only  other  increase  is  in  the  dipole  moment  of  an  average  factor  of  1.25, 
while  the  quadrupole  moment  contribution  is  halved.  The  monopole,  however,  is  26  times 
larger  than  the  dipole,  which  is  the  next  highest  contributing  moment.  The  fourth  through 
the  sixth  moments  are  unaffected  by  the  increase  in  absorption. 

The  presence  of  an  absorption  contrast  enables  the  first  moment  to  be  the  main 
contributor  to  the  magnitude  perturbations.  Higher  levels  of  absorption  contrast  do  not 
increase  the  other  moments  to  be  significant  contributors  in  changing  the  magnitude  of 
the  perturbed  wave. 

Comparing  the  phase  moment  graphs  in  Figures  5.18  and  5.19  shows  that  an  increase 
from  100%  to  200%  contrast  in  absorption  contrast  causes  a  0.1°  increase  in  the  second 
phase  moment.  The  low  scattering  contrast  enables  a  change  in  absorption  contrast  to 
affect  the  phase.  The  remaining  moments  are  unaffected  by  the  increase  in  absorption 
contrast.  However,  the  presence  of  some  absorption  contrast  causes  the  monopole  to  be 
the  major  contributor  to  the  phase  perturbations. 


5-25 


5-4-7  Summary  of  Moment  Analysis.  In  the  absence  of  any  absorption  contrast, 
the  magnitude  and  phase  perturbations  are  dominated  by  the  second  (dipole)  and  third 
(quadrupole)  moments  of  the  scattered  wave.  Adding  an  absorption  contrast  causes  the 
first  moment,  or  the  monopole,  to  be  the  major  contributor  to  both  the  magnitude  and 
phase  perturbations.  In  the  presence  of  a  high  absorption  contrast,  a  low  scattering  contrast 
lowers  the  magnitude  perturbation  contribution  made  by  the  second  moment,  while  an 
increase  in  the  scattering  contrast  causes  a  corresponding  increase  in  the  second,  as  well 
as  the  third  through  sixth  moments.  The  monopole  moment  remains  largely  unaffected  by 
the  presence  of  any  scattering  contrast. 

Presence  of  an  absorption  contrast  affects  the  phase  moments  only  if  there  is  a  low 
scattering  contrast.  In  that  situation,  the  monopole  moment  dominates,  and  the  remaining 
moments  monotonically  decrease  in  contribution  order.  If  scattering  contrast  dominates 
the  phase  perturbations,  the  monopole  moment  falls  to  sixth  rank,  and  the  second  through 
the  fifth  moments  are  the  major  contributors. 
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Figure  5.14  The  first  six  moment  (a)  amplitude  and  (b)  phase  contributions  of  the 
scattered  wave  to  the  incident  wave  are  plotted  for  a  pure  absorber  with 
low  contrast,  along  the  z-axis.  For  both  figures,  fia,out  =  0.05  cm-1, 
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Figure  5.18  The  first  six  moment  (a)  amplitude  and  (b)  phase  contributions  of  the  scat¬ 
tered  wave  to  the  incident  wave  are  plotted  for  high  absorption  and  low 
scattering  contrasts,  along  the  i-axis.  For  both  figures,  na,out  =  0.05  cm-1, 
fJ’s.n-t  =  10  cm-1,  Un  in.  =  0.15  cm-1.  u!.  =  11  cm-1. 


Figure  5.19  The  first  six  moment  (a)  amplitude  and  (b)  phase  contributions  of  the  scat¬ 
tered  wave  to  the  incident  wave  are  plotted  along  the  £-axis.  For  both 
figures,  fia  out  =  0.05  cm-1,  fi'  out  =  10  cm-1,  fia<in  =  0.10  cm-1, 


Figure  5.20  The  first  six  moment  (a)  amplitude  and  (b)  phase  contributions  of  the  scat¬ 
tered  wave  to  the  incident  wave  are  plotted  along  the  i-axis.  For  both 
figures,  fiaout  =  0.05  cm-1,  (i'out  =  10  cm-1,  /x0,m  =  0.10  cm-1, 


5.5  Fourier  Optics  Model  Validity 

Though  the  transfer  functions  of  the  analytic  approach  (Chapter  III)  and  the  Fourier 
optics  wave  propagation  method  (Section  2.5)  cannot  be  compared  directly,  their  models 
of  the  effects  on  the  incident  (or  homogeneous)  wave  due  to  the  presence  of  a  spherical 
inhomogeneity  can  be  observed  and  compared.  The  models  are  simulated  according  to 
the  codes  developed  in  Section  4.8  and  Section  4.9,  respectively.  The  mean  square  error 
(MSE)  is  calculated  in  an  x-y  slice  for  every  0.125  cm  increment  along  the  2-axis  per 
Equation  (4.1).  The  coordinates  corresponding  to  the  inside  of  the  inhomogeneity  are  not 
included  in  the  MSE  calculations  since  the  analytic  transfer  function  is  not  valid  in  that 
region.  The  resulting  MSE  of  the  magnitude  and  phase  of  the  two  approaches  are  plotted 
for  various  absorption  and  scattering  parameters  relative  to  the  background  medium. 

The  Fourier  optics  model  was  developed  under  the  assumption  that  the  wave  prop¬ 
agation  is  through  a  homogeneous,  source-free  medium  [32].  The  inhomogeneity  disrupts 
the  homogeneous  quality  and  subsequently  acts  as  a  secondary  source  within  the  medium. 
Thus  this  model  is  only  valid  in  the  regions  not  containing  the  inhomogeneity  where  the 
homogeneous  and  source-free  conditions  hold.  If  the  inhomogeneity  is  modeled  as  a  thin 
lens  in  the  plane  that  contains  its  center,  then  the  Fourier  optics  model  for  the  incident 
wave  is  valid  up  to  that  plane.  As  the  incident  wave  is  transmitted  through  the  plane 
containing  the  thin  lens  model  of  the  sphere  (Section  2.5.2),  the  incident  wave  is  distorted. 
The  resulting  perturbed  wave  is  then  propagated  through  the  medium  using  the  Fourier 
optics  model  since  the  remaining  medium  is  homogeneous  and  source-free. 

As  seen  in  Figures  5.21  through  5.28,  the  Fourier  optics  model  is  not  without  error. 
The  perturbative  effects  from  the  inhomogeneity  can  only  be  included  as  the  incident  wave 
is  propagated  through  the  plane  containing  the  inhomogeneity,  unlike  the  analytic  transfer 
model  which  can  account  for  perturbations  caused  by  the  inhomogeneity  throughout  the 
entire  medium.  In  a  sense,  there  is  no  history  in  this  forward  propagation  model  as  there 
is  in  the  analytic  transfer  function  model.  Consequently,  the  degree  of  error  between  the 
two  models  increases  as  the  wave  propagates  from  the  point  source  towards  the  plane 
containing  the  center  of  the  inhomogeneity. 


5-35 


Modeling  the  spherical  inhomogeneity  as  a  thin  lens  that  has  straight  ray  propagation 
through  it  rather  than  some  dissipative  behavior,  causes  error  in  both  the  phase  and  the 
magnitude.  Absorption  in  the  sphere  increases  the  amount  of  attenuation  of  the  DPDW 
traveling  through  the  sphere,  while  a  more  scattering  medium  in  the  inhomogeneity  passes 
the  rays  unrefracted.  The  ray  optics  model  does  not  account  for  the  diffraction  behavior. 
As  the  distorted  wave  propagates  away  from  the  plane  containing  the  inhomogeneity  center, 
the  magnitude  is  attenuated,  and  the  MSE  between  the  two  models  will  decrease.  However, 
the  error  introduced  in  the  phase  does  not  attenuate  and  continues  to  propagate  with  the 
perturbed  wave. 

Figures  5.21  and  5.22  depict  the  MSE  for  cases  involving  a  pure  absorber  and  a  pure 
scatterer,  respectively.  The  error  in  both  figures  for  the  homogeneous  wave  up  to  the  plane 
containing  the  center  of  the  inhomogeneity  (z=0)  is  zero  when  there  is  no  contrast  present 
in  either  the  absorption  or  the  scattering  properties.  As  the  amount  of  contrast  increases 
the  amount  of  error  in  both  the  magnitude  and  phase  increases.  The  relative  error  in 
the  phase  is  affected  more  by  an  increase  in  absorption  contrast  than  by  an  increase  in 
scattering  contrast.  Beyond  the  plane  containing  the  inhomogeneity  center,  the  magnitude 
exponentially  drops  off  as  expected.  However,  the  phase  error  caused  by  modeling  the 
inhomogeneity  in  a  plane  is  consistently  carried  through  the  remaining  medium.  A  change 
in  the  absorption  or  scattering  contrast  does  not  significantly  affect  the  phase  error. 

If  a  fixed  absorption  contrast  is  introduced  into  a  pure  scatterer  case,  the  relative 
MSE  for  both  the  magnitude  and  the  phase  increases  in  the  homogeneous  wave  from  the 
source  plane  up  to  the  plane  containing  the  inhomogeneity  center.  The  magnitude  and 
phase  MSE  in  the  perturbed  wave  as  it  propagates  through  the  remaining  medium  is  largely 
unaffected  beyond  the  boundary  of  the  inhomogeneity.  Variation  in  the  scattering  contrast 
does  not  significantly  change  either  the  magnitude  error  or  the  phase  error  in  the  entire 
system.  Figures  5.23,  5.24,  and  5.25  show  the  relative  MSE  for  a  100%,  200%,  and  300% 
fixed  absorption  contrast. 

If  a  fixed  amount  of  scattering  contrast  is  added  to  a  pure  absorber  case,  the  magni¬ 
tude  and  phase  MSE  from  the  source  plane  up  to  the  plane  containing  the  inhomogeneity 
center  is  not  affected  except  when  no  absorption  contrast  is  present.  As  additional  absorp- 
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tion  contrast  is  introduced  into  the  system  as  seen  in  Figures  5.26,  5.27,  and  5.28,  the  MSE 
in  the  homogeneous  wave  propagation  also  increases.  Phase  error  through  the  remaining 
medium,  and  magnitude  error  beyond  the  inhomogeneity  boundary  are  not  affected  by 
variation  in  absorption  contrast. 

These  MSE  results  demonstrate  that  the  Fourier  optics  model  is  valid  for  the  propa¬ 
gation  of  a  spherical  wave  caused  by  a  point  source  in  a  system  containing  an  inhomogeneity 
with  no  optical  contrast  with  the  background  medium.  If  the  presence  of  an  inhomogeneity 
causes  a  contrast,  either  in  absorption  or  scattering,  MSE  in  the  magnitude  of  the  model 
increases  as  the  wave  approaches  the  center  location  of  the  inhomogeneity  and  decreases 
exponentially  as  the  wave  propagates  further  away  from  the  center.  Additional  absorption 
or  scattering  contrast  causes  a  corresponding  error  in  the  magnitude  MSE.  Although  the 
MSE  in  the  phase  of  the  model  also  increases  up  to  the  plane  containing  the  center  of 
the  inhomogeneity,  the  phase  error  caused  by  modeling  the  inhomogeneity  as  a  thin  lens 
remains  in  the  wave  as  it  continues  to  propagate  through  the  medium.  The  MSE  of  the 
phase  is  not  affected  by  additional  contrast  in  either  absorption  or  scattering.  Thus,  the 
Fourier  optics  model  is  a  valid  approximation  to  the  analytic  transfer  function  when  not 
near  the  inhomogeneity,  and  if  the  inhomogeneity  has  a  little  or  no  absorptive  contrast 
with  the  surrounding  turbid  medium.  The  error  incurred  in  the  phase  is  mainly  due  to  the 
model  of  the  inhomogeneity  as  a  spherical  lens,  and  not  an  anomaly  in  the  Fourier  optics 
propagation  model. 
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phases  (degrees)  of  the  total  waves  for  a  pure  absorber  along  the  2-axis 


calculated  analytically  and  by  means  of  the  Fourier  optics  approximations. 


The  object  absorption  coefficient  is  varied  as:  /x0i,n  =  0.05  cm  1  (o), 

Va,in  =  0.10  cm-1  (o),  /zQj in  —  0.15  cm-1  (*),  and  /xa)in  =  0.20  cm-1 
(+)•  For  both  figures,  /j,a  out  =  0.05  cm-1. 


Figure  5.22 


(a) 


(b) 

These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b)  phases 
(degrees)  of  the  total  waves  for  a  pure  scatterer  along  the  2- axis  calculated 
analytically  and  by  means  of  the  Fourier  optics  approximations.  The  object 
scattering  coefficient  is  varied  as:  (isjn  =  10  cm-1  (o),  fis,in  =  11  cm-1 
(°),  Ms, in  =  15  cm-1  (*),  and  Ms, in  =  20  cm-1  (+).  For  both  figures, 
M! 's,out  =  1°  cm_1- 
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Figure  5.23  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b)  phases 
(degrees)  of  the  total  waves  for  an  inhomogeneity  with  low  contrast  in  ab¬ 
sorption  along  the  2-axis  calculated  analytically  and  by  means  of  the  Fourier 
optics  model.  The  object  scattering  coefficient  is  varied  as:  [is,in  =  10  cm-1 
(o),  Ms, in  =  11  cm-1  (o),  nS)in  =  15  cm-1  (*),  and  Ms, in  =  20  cm-1 
(+).  For  both  figures,  n'ain  =  0.10  cm-1,  p a,out  =  0.05  cm-1,  and 

Ms, out  =  10  cm-1. 
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Figure  5.24  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b)  phases 
(degrees)  of  the  total  waves  for  an  inhomogeneity  with  high  contrast  in 
absorption  along  the  z-axis  calculated  analytically  and  by  means  of  Fourier 
optics  model.  The  object  scattering  coefficient  is  varied  as:  ns,in  =  10  cm-1 
(o),  Us, in  =  11  cm-1  (o),  fis>in  =  15  cm-1  (*),  and  fiStin  =  20  cm-1 
(+)•  For  both  figures,  n'ain  =  0.15  cm-1,  na,out  =  0.05  cm-1,  and 

V'sfiut  =  1°  cm-1. 
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Figure  5.25  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b)  phases 
(degrees)  of  the  total  waves  for  an  inhomogeneity  with  very  high  contrast  in 
absorption  along  the  i-axis  calculated  analytically  and  by  means  of  Fourier 
optics  model.  The  object  scattering  coefficient  is  varied  as:  jis,in  =  10  cm-1 
(o),  Ms, in  =  11  cm-1  (o),  nSjin  =  15  cm-1  (*),  and  Ms, in  =  20  cm-1 
(+).  For  both  figures,  \i!ain  =  0.20  cm-1,  Mo, out  =  0.05  cm-1,  and 
Ms, out  =  1°  cm_1- 
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Figure  5.26  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b) 
phases  (degrees)  of  the  total  waves  for  an  inhomogeneity  with  low  contrast 
in  scattering  along  the  5-axis.  The  object  absorption  coefficient  is  varied 
as:  Ha, in  =  0-05  cm-1  (o),  Ha, in  =  0.10  cm-1  (o),  Ha, in  =  0.15  cm-1 
(*),  and  Ha, in  =  0.20  cm-1  (+).  For  both  figures,  H's,in  =  11  cm_1> 
Ha, out  =  0.05  cm-1,  and  Hs,out  =  10  cm-1. 
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Figure  5.27  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b) 
phases  (degrees)  of  the  total  waves  for  an  inhomogeneity  with  high  contrast 
in  scattering  along  the  z-axis.  The  object  absorption  coefficient  is  varied 
as:  Ha, in  =  0.05  cm-1  (o),  Ha, in  =  0.10  cm-1  (o),  Ha, in  =  0.15  cm-1 
(*),  and  Ha, in  =  0.20  cm-1  (+).  For  both  figures,  /z'in  =  15  cm-1, 
Ha, out  =  0.05  cm-1,  and  H's,out  =  10  cm-1. 
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Figure  5.28  These  plots  show  the  mean  square  error  in  the  (a)  magnitudes  and  (b)  phases 
(degrees)  of  the  total  waves  for  an  inhomogeneity  with  very  high  contrast 
in  scattering  along  the  2-axis.  The  object  absorption  coefficient  is  varied 
as:  Ha, in  =  0.05  cm-1  (o),  Ha, in  =  0.10  cm-1  (o),  Ha, in  =  0.15  cm-1 
(*),  and  Ha, in  =  0.20  cm-1  (+).  For  both  figures,  =  20  cm-1, 

Ha, out  =  0.05  cm-1,  and  Hs,out  =  10  cm-1. 
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5.6  Summary 

The  amount  of  perturbation  in  the  total  wave  due  to  the  inhomogeneity  depends  on 
the  contrast  between  the  optical  properties  inside  and  outside  the  object.  A  sensitivity 
analysis  for  various  combinations  of  contrasts  showed  that  an  object  with  a  higher  absorp¬ 
tion  coefficient  relative  to  the  background  medium  causes  a  decrease  in  the  amplitude  of 
the  perturbed  wave  but  does  not  significantly  affect  the  phase.  An  inhomogeneity  with  a 
higher  scattering  parameter  decreases  the  amplitude  of  the  DPDW  but  causes  a  positive 
change  in  the  phase.  A  combination  of  contrasts  combine  their  effects,  meaning  that  the 
reduction  of  the  amplitude  due  to  an  increase  in  absorption  is  in  addition  to  the  reduction 
caused  by  a  contrast  in  scattering  properties.  However,  the  absorption  contrast  does  not 
contribute  significantly  to  a  distortion  in  the  phase,  particularly  in  the  highly  scattering 
cases.  Ultimately,  purely  absorbing  and  scattering  objects  can  be  distinguished  due  to 
the  type  of  change  detected  in  the  DPDW.  An  amplitude  change  indicates  the  presence  of 
an  absorber,  while  a  phase  change  distinguishes  a  scatterer.  However,  when  changes  are 
detected  in  both  amplitude  and  phase,  absolute  characterization  is  not  possible  without 
a  ■priori  knowledge  of  the  object’s  optical  properties. 

The  level  of  perturbation  is  dependent  on  the  location  of  the  detector  plane  to  the 
inhomogeneity.  In  the  simulations  conducted  in  this  research,  the  contrast  causes  the  maxi¬ 
mum  change  in  the  amplitude  and  phase  to  occur  along  the  boundary  of  the  inhomogeneity 
in  the  positive  traveling  direction  of  the  wave.  These  maximum  positions  can  be  used  to 
determine  the  depth  of  the  object.  As  the  diameter  of  the  inhomogeneity  decreases,  the 
corresponding  level  of  perturbation  may  be  much  less  than  what  can  be  actually  detected. 

The  number  of  moments  required  to  characterize  the  scattered  wave  depends  on 
the  level  and  types  of  contrast  present  in  the  system.  In  the  absence  of  any  absorption 
contrast,  the  magnitude  and  phase  perturbations  are  dominated  by  the  second  (dipole) 
and  third  (quadrupole)  moments  of  the  scattered  wave.  Depending  on  the  location  of  the 
detector  plane,  the  fourth  and  fifth  moments  may  contribute  more  than  the  monopole.  If 
an  absorption  contrast  is  introduced  into  the  system,  the  first  moment,  or  the  monopole, 
becomes  the  major  contributor  to  both  the  magnitude  and  phase  perturbations. 
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The  valid  regions  of  the  Fourier  optics  model  of  a  spherical  inhomogeneity  in  an 
otherwise  piece-wise  homogeneous  medium  were  determined  by  comparing  the  total  wave 
perturbations  to  those  calculated  by  the  analytic  transfer  function  model.  The  homo¬ 
geneous  wave  is  best  modeled  by  the  Fourier  optics  model  when  there  is  no  absorption 
contrast  present  in  the  system.  When  there  is  an  absorption  or  a  scattering  contrast 
present,  the  approximate  error  in  the  homogeneous  wave  model  increases  up  to  the  plane 
containing  the  center  of  the  inhomogeneity.  From  that  plane  to  the  boundary  of  the  sys¬ 
tem,  the  magnitude  error  exponentially  decreases,  making  the  wave  propagation  model  a 
more  accurate  model  as  the  detection  plane  moves  away  from  the  inhomogeneity.  However, 
since  the  ray  optics  model  requires  the  spherical  inhomogeneity  to  be  modeled  in  a  plane, 
a  phase  error  is  incurred.  The  wave  propagation  transfer  function  continues  to  propagate 
this  error,  but  does  not  appear  to  incur  any  additional  phase  error.  Effectively,  the  Fourier 
optics  model  is  a  relatively  good  approximation  to  the  analytic  transfer  function  method 
when  not  near  the  inhomogeneity  and  if  the  inhomogeneity  has  a  strong  scattering  and 
low  absorption  behavior.  Otherwise,  the  model  of  the  spherical  inhomogeneity  as  a  thin 
lens  is  not  an  accurate  one. 
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VI.  Summary 


6.1  Conclusions 

Diffusive  light  can  be  used  to  detect  and  localize  optical  inhomogeneities  embedded 
in  thick,  turbid  media  such  as  human  breast  tissue.  In  this  thesis,  a  transfer  function  was 
derived  from  an  exact  analytic  solution  which  can  completely  characterize  the  perturbations 
in  the  forward  propagation  phenomena  caused  by  a  spherical  inhomogeneity  in  a  multiple¬ 
scattering  medium. 

The  diffusion  approximation  to  the  linear  transport  theory  yields  a  wave  solution 
when  the  intensity  of  a  point  source  of  light  is  sinusoidally  modulated.  This  wave  solu¬ 
tion  is  a  highly  damped,  spherically  outgoing  scalar  intensity  wave  which  has  a  complex 
wavenumber.  The  wave  is  referred  to  as  a  diffuse  photon  density  wave  (DPDW)  and 
exhibits  classic  wave  behavior. 

Several  models  have  been  presented  to  describe  the  forward  propagation  of  DPDWs 
through  an  infinite,  homogeneous,  turbid  medium  which  contains  an  embedded  spherical 
inhomogeneity.  Based  on  the  knowledge  of  scattering  of  diffuse  photon  density  waves, 
an  analytic  solution  is  derived  from  the  diffusion  approximation.  This  solution  is  an  ex¬ 
act  model  for  the  perturbations  in  the  DPDW  caused  by  the  inhomogeneity.  Another 
model  uses  Fourier  optics  theory  to  propagate  the  DPDW  through  the  homogeneous  tur¬ 
bid  medium.  This  model  assumes  ray  optic  transmission  through  the  inhomogeneity,  and 
thus  the  model  represents  the  inhomogeneity  as  a  spherical  lens  in  a  plane.  The  DPDW  is 
perturbed  by  the  lens  model  and  is  then  propagated  through  the  remaining  homogeneous 
medium  by  a  simple  plane  wave  transfer  function.  This  lens  model  is  an  approximation 
to  the  perturbations  in  the  DPDW  since  the  ray  optic  assumption  is  not  valid  when  the 
inhomogeneity  exhibits  some  absorptive  property. 

To  improve  the  Fourier  optics  model,  a  transfer  function  was  derived  to  model  ex¬ 
actly  the  perturbations  caused  by  a  spherical  inhomogeneity  in  a  plane  wave  structure. 
This  function  completely  characterizes  the  sensitivity  of  the  system  to  detect  and  local¬ 
ize  in  three-dimensions  inhomogeneities  of  varying  optical  contrasts  with  the  surrounding 
medium.  Through  simulations,  the  sensitivity  analysis  showed  that  a  presence  of  low 
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absorptive  and  scattering  optical  contrasts  can  combine  to  perturb  the  incident  DPDW 
to  a  detectable  level.  The  location  of  the  perturbation  can  then  be  used  to  isolate  the 
inhomogeneity. 

6.2  Recommendations  for  Future  Research 

The  sensitivity  analysis  results  demonstrated  the  remarkable  use  of  perturbative 
behavior  due  to  optical  contrast  as  a  means  to  not  only  detect  but  locate  in  three  dimensions 
a  spherical  inhomogeneity.  However,  these  simulation  results  are  not  complete  without  a 
rigorous  noise  analysis.  The  limitations  of  detectors  are  ultimately  dependent  on  the 
signal-to-noise  ratios  (SNR)  in  the  system.  Previous  SNR  studies  have  been  conducted 
in  systems  that  required  movement  in  both  the  source  and  detector.  The  main  source  of 
error  was  found  to  be  the  source-detector  alignment  [13].  With  the  advent  of  the  exact 
transfer  function  expression  in  the  Fourier  optics  approach  (which  does  not  require  source- 
detector  movement),  the  positional  error  is  removed  from  the  system.  Consequently,  the 
Fourier  optics  approach  has  the  potential  to  achieve  higher  resolution  in  both  detection 
and  localization. 

Perturbations  in  the  DPDW  have  been  shown  to  be  dependent  on  the  level  of  optical 
contrast  present  as  well  as  the  inhomogeneity  size.  Another  parameter  that  has  a  pertur¬ 
bative  effect  and  should  be  investigated  is  the  source  modulation  frequency,  particularly 
in  the  presence  of  low  optical  contrasts.  The  analysis  of  the  frequency  effects  may  also  be 
tied  to  the  SNR  analysis  mentioned  above. 

In  either  the  SNR  or  the  frequency  analysis,  the  conclusions  drawn  from  simulations 
should  be  verified  against  experimental  data.  Initial  experiments  were  conducted  by  Liu 
to  verify  the  theoretical  simulations  [38].  However,  the  SNR  and  frequency  analysis  can 
be  used  to  improve  the  experimental  set-up  to  achieve  higher  accuracy. 

Optical  characterization  continues  to  be  an  elusive  process  [13, 18].  Once  a  tumor 
is  detected,  the  type  of  tumor  becomes  the  important  factor.  Knowledge  of  optical  con¬ 
trasts  may  be  used  to  characterize  detected  tumors.  Low  optical  contrast  was  selected 
as  the  concentration  in  this  thesis  in  anticipation  of  application  to  a  diffraction  tomo- 
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graphic  model  that  is  currently  under  development  [39].  The  analytic  transfer  function 
developed  in  Chapter  III  can  be  directly  applied  to  the  diffraction  tomographic  model  for 
qualitative  reconstruction  of  optical  property  imaging.  Both  theoretical  validation  as  well 
as  experimental  verification  of  the  diffraction  model  using  the  transfer  function  should  be 
pursued. 

6.3  Closing 

The  focus  of  this  research  has  been  on  using  DPDW  imaging  for  low  resolution  breast 
tumor  screening.  This  approach  to  in  vivo  imaging  is  particularly  promising  since  the 
procedure  is  inexpensive  and  has  no  known  adverse  side  effects,  the  noise  associated  with 
movement  of  the  testing  apparatus  is  negligible,  and  the  measurements  can  be  evaluated 
in  a  real-time  environment.  Currently,  industrial  research  groups  are  performing  clinical 
evaluations  of  optical  mammography  systems  [24].  Hopefully  the  results  of  this  thesis  as 
well  as  other  on-going  research  continues  to  generate  interest  and  support  in  the  continued 
development  in  this  and  other  clinical  applications. 
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Appendix  A.  Input  Parameter  Code  Listing  to  the  Analytic  Solution 


This  Appendix  contains  the  source  code  list  for  the  input  parameters  to  the  Analytic 
Solution  code. 


%  Capt  Debbie  Lasocki 

% 

%  file:  M aster. m 

% 

%  Purpose:  Input  parameters 

% 


file_name  =  'Master'; 

%  Constants 

Sac  =  1;  %  Source  modulation  amplitude 

a  =  0.5;  %  Radius  of  inhomogeneity  (cm) 

f  =  20E6;  %  Modulation  fireq.  (Hz) 

numdeg  =  20;  %  Truncate  series  after  this  number  of  terms 

cutoff  =  0.5;  %  Cutoff  frequency  in  filtering  frequency  domain  data 

%  Geometry  of  Source  and  location  of  Detector  plane  (cm) 

%  relative  to  the  center  of  the  inhomogeneity  (0,0,0) 

rsx  =  0;  %  Location  of  the  source  along  x-axis 

rsy  =  0;  %  Location  of  the  source  along  y-axis 

rsz  =  —2;  %  Location  of  the  source  along  z-axis 

zdet  =  3;  %  Location  of  detector  plane 


%  Define  the  region  of  interest  and  the  increments  of  resolution  in 
%  units  of  centimeters.  Do  not  include  origin  in  the  sampling. 


inc  =  0.1; 
idim  =  5.0; 
zinc  =  0.125; 


x  — 

y  = 

z  = 


— idim/2:inc:idim/2]; 

— idim/2:inc:idim/2]; 
rsz+zinc:zinc:— zinc  zinc:zinc:zdet]; 


xplot  =  ceil(length(x)/2);  %  y-plane  at  which  to  view  x-z  slice  (pixel) 

yplot  =  ceil(length(y)/2);  %  y-plane  at  which  to  view  x-z  slice  (pixel) 

zplot  =  length(z);  %  z-plane  at  which  to  view  x-y  slice  (pixel) 

%  Parameters  outside  of  the  inhomogeneity 

Muaout  =  0.05;  %  absorption  coeff  (cm~-l) 

Musout  =  10;  %  red.  scattering  coeff  (cm‘-l) 

Nout  =  1.333;  %  index  of  refraction 

Vout  =  2.98E10/Nout;  %  speed  of  light  without  medium 

Dout  =  Vout  /  (3*(Musout+Muaout));  %  Diffusion  coeff 

Kout  =  sqrt((— Vout*Muaout  +  2*pi*f*j)/Dout);  %  wavenumber 


%  Parameters  inside  of  the  inhomogeneity 
Muain  =  0.15; 

Musin  =  12; 

Nin  =  1.333; 

Vin  =  2.98E10/Nin; 

Din  =  Vin  /  (3*(Musin+Muain)); 

Kin  =  sqrt((— Vin*Muain  +  2*pi*f*j)/Di: 


%  absorption  coeff  (cm~-l) 

%  red.  scattering  coeff  (cm~-l) 
%  index  of  refraction 
%  speed  of  light  within  medium 
%  Diffusion  coeff 
l);  %  wavenumber 
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Appendix  B.  Analytic  Solution  Code  Listing 

This  Appendix  contains  the  source  code  list  for  the  Analytic  Solution. 


function  [Uinc,Uscatt,Us,Uout,UoutH,H]  =  trx(numdeg,  . . . 

Sac,  . . . 
a,  . . . 

rsx,  rsy,  rsz,  . . . 
x,  y,  z,  ... 

Dout,  . . . 

Vout,  . . . 

Kout,  . . . 

Din,  . . . 

Kin  . . . 

); 

%  Gapt  Debbie  Lasocki 

% 

%  file:  trx.m 

% 

%  Purpose:  To  examine  the  Scattered  wave  as  derived  from  an  analytic  solution 

% 

%  It  is  assumed  that  the  incident  wave  is  generated  by  a  point  source  and 
%  the  detector  is  in  a  plane  on  the  opposite  side  of  the  inhomogeneity  than 
%  the  source.  Further,  the  inhomogeneity  is  set  at  the  origin  of  this 
%  coordinate  system,  and  all  distances  are  measured  relative  to  that  point. 

%  This  can  be  set  up  to  calculate  the  wave  at  various  z-planes,  but  the 
%  planes  of  interest  here  are  already  calculated  meaning  that  at  the 
%  inhomogeneity  center  and  at  the  detector  plane,  the  wave  values  are  already 
%  computed. 

%  Calculate  dimensions 

lx  =  length(x); 
ly  =  length(y); 

Iz  =  length(z); 

%  Determine  the  Incident  Wave 

%  The  incident  wave  is  generated  by  a  point  source  in  the  homogeneous  medium, 
%  but  outside  the  region  of  interest  so  singularities  are  avoided. 

%  Pre-calculate  the  x-y-z  distances  in  spherical  coordinates  about  the 
%  source  coordinate 
[xp,yp,zp]  =  ndgrid(x,y,z); 

R  =  sqr  t  ( (xp — rsx) .  ~  2  +  (yp-rsy)."2  +  (zp-rsz).‘2); 

Uinc  =  Vout*Sac./(4*pi*Dout.*R)  .*  exp(j*Kout.*R); 

%  A  (1,0)  coefficients  defined 
for  l=0:numdeg— 1 

terml(l+l)  =  — j*Sac*Kout*sbesselh(l,l,Kout*abs(rsz)). . . 
*sqrt((2*l+l)/(2*pi))*conj(ymn('e '  ,0,l,pi,0)); 

term2(l+l)  =  Dout. *Kout.*a.*Dsbesselj(l, Kout. *a).*sbesselj(l, Kin. *a). . . 

—  Din.  *Kin.*a.*sbesselj(l, Kout.  *a).*Dsbesselj(l, Kin.  *a); 

term3(l+l)  =  Dout. *Kout.*a.*Dsbesselh(l,l, Kout. *a).*sbesselj(l. Kin. *a). .. 

—  Din.  *Kin.*a.*sbesselh(l,l,  Kout.  *a).*DsbesseIj(l,Kin.*a); 

end; 

Aim  =  terml  .*  term2./term3; 

%  Define  scattered  wave  in  the  region  of  interest. 

%  Pre-calculate  the  x-y-z  distances  in  spherical  coordinates  about  the 
%  detector  plane  wrt  to  the  origin  of  the  inhomogeneity 

R  =  sqrt(xp.*2  +  yp.*2  +  zp.“2); 
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for  l=0:numdeg— 1; 
for  xp  =  l:lx 
for  yp  =  l:ly 
for  zp  =  l:lz 


%  For  each  moment 
%  For  each  x-point 
%  For  each  y-point 
%  For  each  z-point 


if  R(xp,yp,zp)  '=  0 

Theta  =  acos(z(zp)./R(xp,yp,zp)); 

else 

Theta  —  0; 

end; 

if  x(xp)  "=  0 

Phi  =  atan(y(yp)./x(xp)); 

else 

Phi  =  pi/2; 
end; 

Uscatt(xp,yp,zp,l+1)  =  Alm(l+1)  ... 

*  sbesselh(l,l,Kout*R(xp,yp,zp))  . . . 

*  sqrt((2*l+l)/(2*pi))*ymn( 1  e 1 ,0,l,Theta,Phi); 

end; 

end; 

end; 

end; 


Us  =  sum(Uscatt,4);  %  Sum  over  the  moments 


%  Output  wave  is  just  the  superposition  of  the  incident  and  the 
%  scattered  wave 


Uout  =  Uinc  +  Us; 
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Appendix  C.  Transfer  Function  Code  Listing  to  the  Analytic  Solution 

This  Appendix  contains  the  source  code  list  for  the  transfer  function  to  the  Analytic 
Solution  code.  The  first  code  is  the  main  program  which  calls  a  transfer  function  calculation 
which  is  listed  in  the  second  code. 


function  [Uinc,Uscatt,H,UoutH]  =  fth(numdeg,  . . . 

Sac,  ... 
a,  . . . 

rsx,  rsy,  rsz,  . . . 
x,  y,  z,  ... 

Dout,  . . . 

Vout,  ... 

Kout,  ... 

Din,  . . . 

Kin  ...  10 


%  Capt  Debbie  Lasocki 

% 

%  file:  fth.m 

% 

%  Purpose:  To  examine  the  Scattered  wave  as 
%  determined  by  a  transfer  function 

%  derived  from  an  analytic  solution 

% 

%  It  is  assumed  that  the  incident  wave  is  generated  by  a  point  source  and 
%  the  detector  is  in  a  plane  on  the  opposite  side  of  the  inhomogeneity  than 
%  the  source.  Further,  the  inhomogeneity  is  set  at  the  origin  of  this 
%  coordinate  system,  and  all  distances  are  measured  relative  to  that  point. 

%  This  can  be  set  up  to  calculate  the  wave  at  various  z-planes,  but  the 
%  planes  of  interest  here  are  already  calculated  meaning  that  at  the 
%  inhomogeneity  center  and  at  the  detector  plane,  the  wave  values  are  already 
%  computed. 

%  Calculate  dimensions 

lx  =  length(x); 
ly  =  length(y); 
lz  =  length(z); 

xcent  =  ceil(lx/2);  %  center  of  x-plane  (pixel) 

ycent  =  ceil(ly/2);  %  center  of  y-plane  (pixel) 

%  Determine  the  Incident  Wave 

%  The  incident  wave  is  generated  by  a  point  source  in  the  homogeneous  medium, 
%  but  outside  the  region  of  interest  so  singularities  are  avoided. 

%  Pre-calculate  the  x-y-z  distances  in  spherical  coordinates 
%  about  the  source  coordinate 
[xp,yp,zp]  =  ndgrid(x,y,z); 

R  =  sqrt((xp— rsx). "2  +  (yp— rsy). ‘2  +  (zp— rsz). "2); 

Uinc  =  Vout*Sac./(4*pi*Dout.*R)  .*  exp(j*Kout.*R); 

%  A  (1,0)  coefficients  defined 
for  l=0:numdeg— 1 

terml(l+l)  =  — j*Sac*Kout*sbesselh(l,l,Kout*abs(rsz)). . . 
*sqrt((2*l+l)/(2*pi))*coiy(ymn(le,,0,l,pi,0)); 

term2(l+l)  =  Dout. *Kout.*a.*Dsbesselj(l, Kout. *a).*sbesselj(l,Kin.*a). .. 

—  Din. *Kin.*a.*sbesselj(l, Kout. *a).*Dsbesselj(l,Kin.*a); 

term3(l+l)  =  Dout. *Kout.*a.*Dsbesselh(l,l, Kout. *a).*sbesselj(l, Kin. *a). .. 
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—  Din.*Kin.*a.*sbesselh(l,l,Kout.*a).*Dsbesselj(l,Kin.*a); 

end; 

Aim  =  terml  .*  term2./term3; 

%  Determine  the  scattered  wave  in  the  region  of  interest 


%  Pre-calculate  the  x-y-z  distances  in  spherical  coordinates  about  the 
%  detector  plane  wrt  to  the  origin  of  the  inhomogeneity 

R  =  sqrt(xp.*2  4-  yp.*2  +  zp."2); 


for  l=0:numdeg— 1; 
for  xp  =  l:lx 
for  yp  =  l:ly 
for  zp  =  l:lz 


%  For  each  moment 
%  For  each  x-point 
%  For  each  y-point 
%  For  each  z-point 


if  R(xp,yp,zp)  '=  0 

Theta  =  acos(z(zp)./R(xp,yp,zp)); 

else 

Theta  =  0; 

end; 

if  x(xp)  *=  0 

Phi  =  atan(y(yp)./x(xp)); 
else 

Phi  =  pi/2; 

end; 

Uscatt(xp,yp,zp,l+1)  =  Alm(l+1)  . . . 

*  sbesselh(l,l,Kout*R(xp,yp,zp))  ... 

*  sqrt((2*l+l)/(2*pi))*ymn(,e' ,0,1, Theta, Phi); 

end; 

end; 

end; 

end; 


Us  =  sum(Uscatt,4);  %  Sum  over  the  moments 

%  Determine  Transfer  function 

%  Zero  pad  up  to  automatically-defined  fftpts  size  to  avoid  wrap-around  and 
%  fftshift  problems  on  an  odd-sized  array  (in  x-y  direction). 


%  The  center  of  the  medium  must  be  shifted  to  the  upper  left-hand  comer  of 
%  the  array  in  order  to  use  the  fft  command  properly  in  MATLAB.  A  second 
%  shift  is  done  to  return  the  array  to  the  proper  ordering  and  the  original 
%  data  size  is  extracted  out  of  the  center  of  the  array. 


fftpts  =  2~ceil(log2(lx)); 
tine  =  zeros(fftpts, fftpts, lz); 
ts  =  zeros(fftpts,  fftpts,  lz); 
fftcent  =  fftpts/2  +1; 
xi  =  — xcent+l:xcent— 1; 
yi  =  — y  cent +l:y  cent— 1; 


%  Assumes  x-axis  is  longest 
%  initialize  temp  array 
%  initialize  temp  array 
%  Center  of  fft  matrix 
%  midrange  array  indices  x-axis 
%  midrange  array  indices  y-axis 


tinc(fftcent+xi,  fftcent+yi,l:lz)  =  Uinc(:,:,:); 
ts(fftcent+xi,  fftcent+yi,l:lz)  =  Us(:,:,:); 


for  n  =  l:lz; 

UincFT(:,:,n)  =  fftshift(ifft2(fftshift(tinc(:,:,n)))); 
UscattFT(:,:,n)  =  fftshift(ifft2(fftshift(ts(:,:,n)))); 

end; 


%  A  matrix  of  the  same  size  of  the  fft  arrays  is  created  in  order  to 
%  correctly  add  1  to  each  pixel. 


unity  =  ones(fftpts, fftpts, lz); 

H  =  unity  +  UscattFT./UincFT; 


%  Determine  Output  wave  from  Transfer  function 
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%  The  shift  is  done  to  correctly  order  the  array  for  fft. 

%  The  second  shift  is  done  at  plotting  the  results  to  extract 

%  desired  original  data  (slice)  from  the  results  of  the  transform. 

for  n=l:lz 

UH(:,:,n)  =flftshift(fft2(flftshift(UincFT(:,:,n).*H(:,:,n)))); 
end; 

UoutH  =  zeros  (lx, ly,lz); 

UoutH  =  UH(fftcent+xi,  fftcent+yi,:); 
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Appendix  D.  Fourier  Optics  Wave  Propagation  Code  Listing 

This  Appendix  contains  the  source  code  list  for  the  Fourier  optics  wave  propagation  model 
[41]  using  Goodman’s  methods.  The  first  code  is  the  main  program  which  calls  a  transfer 
function  calculation  which  is  listed  in  the  second  code. 

function  [homo.Udet]  =  goodman(  Sac,  . . . 

a,  . . . 

rsx,  rsy,  rsz,  . . . 
inc,  . . . 
idim,  . . . 
zinc,  . . . 
x,  y,  z,  ... 

Dout,  . . . 

Vout,  ... 

Kout,  ...  10 

Kin  . . . 

); 

%  Capt  Debbie  Lasocki 

% 

%  file:  goodman.m 

% 

%  Purpose:  To  propagate  a  spherical  wave  through  an  otherwise 
%  piece-wise  homogeneous  media  containing  a  spherical  object 

%  using  plane  wave  propagation  methods  developed  by  Goodman.  20 

%  This  program  was  ported  to  MATLAB  5.0  from  IDL  code  developed  by 
%  Dr.  Chuck  Matson. 

%  This  program  calculates  an  incident  wave  through  a  a  homogeneous  media. 

%  It  also  calculates  the  wave  through  the  homogeneous  media  up  to  the 
%  center  of  an  inhomogeneity  via  Goodman’s  approach.  The  wave  is  then 
%  transmitted  through  the  media  using  conventional  techniques.  This 
%  perturbed  wave  is  then  propagated  to  the  detector  plane  using  Goodman’s 

%  approach  once  again.  30 

%  The  wave  is  centered  on  the  x-y  axis  and  symmetric  around  the  z-axis. 

%  The  resolution  of  propagation  along  each  axis  is  determined  by  the  user. 

%  The  size  of  the  sample,  the  location  of  the  source,  the  location  of  the 
%  inhomogeneity,  and  the  optical  parameters  of  the  inhomogeneity  are  also 
%  specified  by  the  user. 

%  Calculate  dimensions 

ic  =  idim/2;  %  Center  of  sample  (cm)  40 

flr  =  IE— 10;  %  Establish  floor  value 

lx  =  length(x);  %  length  of  x-axis  (pixels) 

ly  =  length(y);  %  length  of  y-axis  (pixels) 

lz  =  length(z);  %  length  of  z-axis  (pixels) 

lzsph.  =  find(z==0);  %  center  of  sphere  along  z-axis  (pixels) 

%  assuming  source  pixel  is  included  in  z 

%  Size  of  FFT  matrix  should  be  set  >  It  pixels  in  sample  50 

fftpts  =  2*ceil(log2(lx));  %  Assumes  x-axis  is  longest 

%  Pre-calculate  the  x-y-z  distances  in  spherical  coordinates  about  the 
%  source  coordinate 
[X  Y  Z]  =  ndgrid(x,y,z); 

XY2  =  (X— rsx).‘2  +  (Y-rsy).*2; 

R  =  sqrt(XY2  +  (Z-rsz).‘2); 

%  Generate  incident  wave  in  the  detection  plane 
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homo  =  Sac*Vout./(4*pi*Dout.*R).*exp(j*Kout.*R); 
size(homo) 

%  Propagate  source  wave  to  the  plane  containing  center  of  inhomogeneity 
%  which  also  corresponds  to  the  origin  of  the  system 
%  This  creates  a  separate  incident  wave  then  the  one  above  (homo). 

%  The  field  is  generated  at  each  z-plane  increment  up  to  the  center 
%  of  the  inhomogeneity  along  the  z-axis 

Uhomo  =  Sac*Vout./(4*pi*Dout.*R(:,:,l:lzsph)).*exp(j*Kout.*R(:,:,l:lzsph)); 

70 

%  Create  the  complex  transmittance  caused  by  the  sphere 
%  All  other  points  outside  the  sphere  are  transmitted  unperturbed 

%  Determine  location  of  center  of  sphere  in  pixels 
icenx  =  ceil(lx/2); 
iceny  =  ceil(ly/2); 

if  rem(icenx,2)  ==  0 
ix  =  icenx; 

else  80 

ix  =  floor(icenx); 
end; 

if  rem(iceny,2)  ==  0 
iy  =  iceny; 
else 

iy  =  floor(iceny); 
end; 

trans  —  ones(lx);  %  Initialize  complex  transmittance  array 

90 

for  m  =  — a/inc:a/inc  %  Index  array  for  diameter  of  sphere  in  pixels 

for  n  =  — a/inc:a/inc 

if  a"2  >=  XY2(ix+m,iy+n)  %  Is  point  inside  sphere? 

thick  =  2.*sqrt(a*2— XY2(ix+m,iy+n));  %  Distance  through  sphere 
amp  =  exp(— thick*imag(Kin— Kout));  %  Adjust  amplitude 

phase  =  thick* (Kin— Kout);  %  Adjust  phase 

trans(ix+m,iy+n)  =  amp*exp(j*phase);  %  Calculate  transmittance 
end; 
end; 

end;  100 

%  Create  perturbed  wave  in  same  plane,  just  after  going  through  equivalent 
%  sphere  transmittance  as  created  above  (so  plane  location  is  at  the  center 
%  of  the  inhomogeneity  still) 

Uper=Uhomo; 

Uper(:,:,lzsph)  =  Uhomo(:,:,lzsph).*trans; 

%  Propagate  inhomogeneous  wave  to  detection  plane  now  using  Goodman’s 

%  Fourier  optics  propagation  approach  110 

%  Zero  pad  perturbed  wave  to  avoid  wrap-around  errors  and  quadrant  swapping 
%  errors  during  Fourier  transform  process 
Uperz  =  zeros(fftpts,fftpts,lzsph); 

UPERZ  =  zeros(fftptSjfftpts,lzsph); 

fftcent  =  fftpts/2  +1;  %  Center  of  fft  matrix 

xi  =  — icenx+l:icenx— 1;  %  midrange  array  indices  x-axis 

yi  =  — iceny+l:iceny— 1;  %  midrange  array  indices  y-axis 

Uperz(fftcent+xi,  fftcent+yi,l:lzsph)  =  Uper(:,:,:); 

120 

for  n  =  l:lzsph 

%  Fourier  transform  perturbed  wave 
UPERZ(:,:,n)  =  fFtshift(fFt2(fFtshift(Uperz(:,:,n)))); 

end; 

%  Initialize  arrays 
H  =  ones(fftpts,  fftpts,  Iz); 


D-2 


130 


UH  =  zeros(fftpts,  fftpts,  lz); 

UH(:,:,l:lzsph)  =  UPERZ(:,:,:); 
uh  =  zeros(fftpts,  fftpts,  lz); 

%  Pre-calculate  transfer  function 
Htemp  =  transfer(fftpts,  zinc,  inc,  Kout,  fir); 

for  n=lzsph+l:lz; 

%  Move  wave  via  transfer  function 
H(:,:,n)  =  H(:,:,n).*Htemp; 

UH(:,:,n)  =  UH(:,:,n-l).*H(:,:,n); 

end; 

140 

for  n  =  l:lz 

%  Inverse  Fourier  transform,  to  get  inhomogeneous  wave 
uh(:,:,n)  =  fftsh!ft(ifft2(fftshift(UH(:,:,n)))); 

end; 

Udet  =  zeros(lx,ly,lz); 

Udet  =  uh(fftcent+xi,fftcent+yi,:);  %  Extract  original  size  of  wave  array 

150 


function  [H]  =  transfer(fftpts,  zprop,  scale,  k,  floorratio); 

%  This  function  creates  a  (Matson)  transfer  function  for  wave  propagation. 

%  fftpts  is  the  size  of  the  array,  zprop  is  the  distance  to  propagate  between 
%  orthogonal  planes,  scale  is  the  distance  per  pixel  in  the  image  plane, 

%  and  k  is  the  wave  number  (complex),  and  the  floorratio  is  the  lower 
%  limit  for  values  to  keep  within  the  precision  of  MATLAB  5.0. 

% 

%  fftpts  :  size  of  array  10 

%  zprop  :  plane  of  interest  distance 
%  scale  :  x-y  increments  in  image  plane 
%  k  :  complex  wavenumber 

%  floorratio  :  bottom  limit  of  transfer  function 

ic  =  fftpts/2  +  1;  %  Center  of  array 

%  Note  that  this  is  one-pixel  short  of  the 
%  center  of  the  fftpts  array 
%  Results  in  H  are  centered  on  the  fftpts 

%  array  as  it  should  be  20 

H  =  zeros  (fftpts);  %  Initialize  array 

fscale  =  2*pi/(fftpts*scale); 

m  =  l:fftpts;  n  =  l:fftpts; 

[M  N]  =  meshgrid(m,n); 

r  =  j*zprop*sqrt(k.'2  —  (fscale*(ic— M))."2  —  (fscale*(ic— N)).*2); 
rmag  =  abs(r); 

rphasor  =  r./rmag;  30 

r  =  rmag.*rphasor; 

H  =  exp(r); 

H  =  H+floorratio  *  max(max(H))*ones(size(H)); 
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